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Introduction 


Early  detection  of  masses  within  the  breast  that  may  transform  into  malignancies  is  known  to  be 
essential  for  positive  treatment  and  outcome.  Currently,  mammography  is  the  clinical  standard 
for  screening  and  provides  useful  but  at  times  ambiguous  information,  which  can  necessitate 
further  invasive  workup  of  benign  lesions.  Recent  research  has  indicated  that  a  family  of 
methods  developed  in  the  field  of  elastography  holds  promise  as  an  alternative  means  of 
interrogating  soft  tissue  structures  by  providing  a  spatial  mapping  of  material  properties  (e.g. 
elasticity)  that  can  be  inspected  for  the  detection  of  lesions  [1,  2],  A  technique  known  as 
‘modality-independent  elastography’  (MIE)  [3-5]  has  been  introduced  with  the  intent  of 
combining  the  intuitive  discrimination  from  manual  detection  with  the  superior  depth  of 
penetration  and  anatomical  detail  typically  given  by  imaging.  The  basic  requirements  for  the 
method  are  two  images  of  the  tissue  in  different  states  of  deformation  (e.g.  compression). 
Elasticity  parameters  are  then  reconstructed  within  the  context  of  an  inverse  problem  that  utilizes 
non-rigid  image  registration  constrained  by  a  biomechanical  model  in  order  to  best  describe  the 
composition  of  the  tissue.  The  final  result  is  a  map  of  the  breast  (or  other  tissue  of  interest)  that 
reflects  material  inhomogeneity,  such  as  in  the  case  of  a  tumor  mass  that  disrupts  the  surrounding 
structure  of  normal  tissue.  Because  MIE  works  on  probing  the  differences  between  images,  it 
can  be  used  to  not  only  work  in  concert  with  more  traditional  screening  techniques  but  also 
address  a  possible  gap  when  those  methods  are  unable  to  directly  discern  tissues  of  interest. 


Body 

As  stated  in  the  original  proposal,  there  are  three  main  aims  of  this  project:  (1)  to  expand  and 
refine  the  current  MIE  technique  to  enhance  its  efficiency  and  capabilities,  (2)  to  perform 
analyses  on  texture  in  input  images  and  quantify  statistical  parameters  capable  of  estimating  and 
evaluating  the  success  of  elastographic  reconstruction,  and  (3)  to  engineer  a  device  that  is 
compatible  with  current  medical  imaging  systems  and  can  produce  compressive  forces 
appropriate  for  phantom  and/or  clinical  setups.  The  relevant  proposed  and  completed  work  is 
listed  below  and  organized  as  closely  as  possible  to  the  two  major  arcs  of  the  Statement  of  Work. 

Task  1(a)  stated:  “Incorporation  of  additional  biomechanical  models  (e.g.  3D  and  2D/3D 
deformation  effects).” 

Initial  work  with  MIE  involved  reconstructions  of  circular  stiff  inclusions  embedded  within  thin 
rectangular  membranes  of  polyurethane  rubber.  The  materials  used  have  essentially 
indistinguishable  colors  but  vary  significantly  in  their  elastic  modulus  values  (a  contrast  ratio  of 
approximately  5.7:1  as  confirmed  by  material  testing  by  an  Enduratec  ELF  3200).  A  black 
permanent  marker  was  used  to  place  a  pattern  of  regularly  spaced  (~1  cm)  grid  lines  across  the 
membrane,  which  was  then  securely  clamped  along  two  opposite  edges  and  subjected  to  a 
uniaxial  tensile  displacement  (~8%  strain)  by  means  of  a  milling  vise.  A  commercial  webcam 
was  rigidly  mounted  above  the  membrane  to  acquire  image  pairs  of  the  pre-  and  post¬ 
deformation  states.  The  intent  of  these  experiments  was  to  have  a  means  of  performing  relatively 
quick  2D  analyses  and  possibly  develop  an  application  for  dermoscopic  evaluation  of  skin 
lesions.  The  model  utilized  the  plane  stress  approximation  [6],  which  is  actually  a  reduction  of 
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the  3D  formulation  of  Cauchy’s  Law.  However,  because  the  skin  has  multiple  layers  in  which  a 
potential  lesion  might  be  present,  it  was  uncertain  whether  the  use  of  planarity  was  suitable  for 
elastographic  work.  In  particular,  the  subcutaneous  layer  of  adipose  tissue  that  serves  to 
essentially  lubricate  the  epidermis  and  dermis  may  have  a  non-trivial  thickness  such  that  in-plane 
loading  is  no  longer  confined  to  a  sufficiently  thin  slab.  Therefore,  a  2D/3D  comparison  was 
performed  to  examine  the  effects  of  model  assumptions  on  deformations  used  in  the 
reconstruction  process.  Two  physical  characteristics  of  the  domain  were  examined:  the  depth  of 
penetration  of  the  stiff  lesion  (simulating  a  melanoma)  and  the  thickness  of  the  subcutaneous 
layer.  Full  3D  deformations  were  solved  over  a  finite  element  mesh,  projected  to  the  surface, 
and  aligned  to  the  dermoscopic  image  before  proceeding  to  the  elastographic  reconstruction. 
Figures  1  and  2  show  the  coupling  of  the  3D  model  with  2D  reconstruction.  The  analysis  of  the 
recovered  elasticity  distributions  indicated  that  deeper  lesions  were  more  accurately 
characterized,  and  the  thickness  of  the  subcutaneous  layer  had  no  significant  impact  on  the 
reconstruction.  Further  details  of  this  work  can  be  found  in  Appendix  D. 
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Figure  1.  MIE  2D/3D  experiment  setup,  (a)  finite  element  mesh  (b)  overlay  of  displacements  induced  by 
hypothetical  dermoscopic  probe  (c)  diagram  of  skin  composition  in  mesh  (gray=lesion)  from  side  view. 
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Figure  2.  MIE  2D/3D  images,  (a)  source  image  obtained  from  a  digital  photograph  of  a  melanoma  (b)  sample 
simulation  reconstruction  obtained  by  experiment  shown  in  Figure  1.  Black  overlays  in  both  frames  denotes  border 
of  lesion. 


Task  1(b)  stated:  “Design  and  implementation  of  computational  code  routines  (e.g. 
modification  of  finite  element  model).”  and  Task  1(c)  stated:  “Enhancement  of  model 
programming  structure  via  numerical  analysis  and  parallelization.” 

In  order  to  accommodate  the  methodology  of  MIE  in  creating  a  Jacobian  matrix  fully  sensitive  to 
the  discretization  of  the  domain,  a  large  number  of  solutions  involving  the  finite  element  model 
and  the  subsequent  imge  deformation  are  required.  With  the  proposed  increase  in  dimensionality 
to  handle  3D  data,  the  implementation  complexity  quickly  increases  beyond  the  capabilities  of 
the  original  MATLAB/FORTRAN/LAPACK  design.  Therefore,  the  Portable  Extensible  Toolkit 
for  Scientific  Computation  (PETSc)  toolkit  [7,8]  was  selected  to  provide  the  necessary 
framework  for  developing  sparse  matrix  system  solvers  and  split  the  Jacobian  formation  process. 
A  separate  C/C++  routine  has  also  been  written  to  perform  a  Gauss-Newton  optimization  and 
interface  with  PETSc  solver  structures.  This  framework  is  flexible  enough  to  accommodate  the 
use  of  2D  models  such  that  unified  data  structures  and  formats  are  now  employed,  and  multi¬ 
resolution  capabilities  present  in  the  previous  2D-only  codebase  are  now  available  in  3D. 

It  was  previously  reported  that  by  using  a  share  of  100  CPUs  from  the  Vanderbilt 
University  Advanced  Computing  Center  for  Research  and  Education  cluster,  speedups  in  the 
overall  computation  of  over  two  orders  of  magnitude  were  achieved.  This  actually  indicated  a 
superlinear  path  over  predictions  from  Ahmdal’s  Law  [9];  however,  recent  observations  indicate 
that  although  significant  gains  are  regularly  made  in  employing  parallel  code,  heterogeneity  of 
machine  power  and  network  latencies  may  actually  hinder  efficiency.  This  is  hypothesized  to  be 
a  consequence  of  the  current  design  of  the  internal  messaging  system  and  may  be  resolved  by 
adjusting  the  current  master-slave  relationship  to  utilize  a  dynamic  SPMD  (single  program, 
multiple  data)  flow  concept  wherein  the  master  attempts  to  concurrently  designate  a  task  to  the 
next  available  processor  while  accepting  the  latest  completed  job.  An  in-house  cluster  of  18 
processors  (2.0  Ghz  Pentium  4  Xeon)  is  now  available  for  performing  additional  reconstructions 
and  further  benchmarking  to  investigate  this  enhancement. 
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Task  1(d)  stated:  “Design  of  texture  analysis  parameter  for  predicting  elastographic 
reconstruction  effectiveness.” 

Earlier  work  had  examined  the  effects  in  2D  of  noise  in  both  the  image  and  boundary  condition 
selection.  The  primary  conclusions  of  that  work  were  that  degradation  in  the  image  of  up  to  10% 
randomized  noise  was  within  the  expected  behavior  of  most  acquisition  systems,  and  that 
boundary  mismappings  of  less  than  0.5  pixel  units  was  generally  tolerated  by  the  reconstruction 
algorithm.  For  further  detail,  please  see  the  full  text  of  this  work  as  listed  in  Appendix  B. 

With  the  recent  move  to  3D  reconstruction  work,  it  has  become  of  increasing  concern  to 
analyze  the  effects  of  boundary  condition  estimation.  Because  similarity  metrics  (e.g.  correlation 
coefficients)  typically  normalize  and/or  blur  the  intensity  distribution  [10],  image  noise  can  be 
tolerated  fairly  well.  However,  any  inaccuracy  in  the  boundary  condition  set  propagates  through 
the  model  and  into  image  deformation,  thereby  having  a  significant  impact  on  the  evaluation  of 
the  objective  function  and  the  success  of  the  reconstruction.  The  current  method  of  selecting 
displacement  (Type  1)  boundary  conditions  requires  manual  interaction  to  guide  or  correct  point 
correspondence  for  every  surface  node.  Assuming  that  an  input  device  (e.g.  a  mouse)  is  needed 
to  identify  the  specific  coordinates,  this  introduces  an  operator-dependent  noise  process  in 
localizing  any  given  point.  A  gold  standard  boundary  condition  set  was  derived  by  simulating 
compression  of  the  breast  with  a  device  described  below  in  the  work  related  to  Task  2. 
Reconstructions  were  performed  using  randomized  vectors  to  disrupt  the  target  surface;  the 
results  of  this  experiment  indicate  that  improper  localization  of  boundary  points  greater  than  or 
equal  to  0.5  units  of  voxel  spacing  can  introduce  significant  error  to  the  reconstruction  process 
and  impair  its  ability  to  characterize  the  underlying  elasticity  distribution.  This  is  a  similar  result 
to  that  found  from  the  prior  work  done  in  2D  [Appendix  B],  It  also  confirms  that  randomizing 
the  vectors  is  a  significant  challenge  to  the  algorithm  because  of  the  introduction  of  highly  non¬ 
physical  deformations  that  cause  backlash  in  the  finite  element  mesh  and  other  numerical 
anomalies. 

The  increase  in  dimensionality  from  2D  to  3D  necessitates  much  higher  discretization 
levels  in  meshing  and  pre-processing  work.  To  effectively  cover  the  surface  of  a  CT  volume 
requires  over  6,000  nodes,  which  would  result  in  an  unreasonable  task  for  manual  boundary 
condition  selection.  Therefore,  an  automated  method  for  determining  surface  correspondence  and 
then  interpolating  those  results  into  displacements  for  boundary  condition  sets  is  needed.  Three 
methods  of  surface  registration  and  point  correspondence  were  considered;  two  are  derived  from 
surface  matching  of  potential  energy  distributions  based  on  the  diffusion  and  Laplace  equations, 
and  the  other  is  a  ffee-form  warping  via  a  thin-plate  spline.  For  complete  details  on  the  methods 
and  buildup  from  previous  work,  please  see  documents  provided  in  Appendix  D  and  Appendix 
F.  ompared  to  the  gold  standard  boundary  condition  set,  the  target  registration  error  was  greatest 
for  the  diffusion  method  and  least  (best)  for  the  thin-plate  spline.  Energy  matching  from  the 
solutions  of  the  diffusion  and  Laplace  equations  yielded  boundary  condition  sets  that  were 
inadequate  for  reconstructing  a  proper  elasticity  contrast.  This  can  be  partly  explained  because 
the  mean  errors  of  those  surface  registration  techniques  (as  compared  to  the  gold  standard  and  in 
an  equivalent  sense  to  the  tested  noise  levels)  are  approximately  3.3  and  1.7  voxel  units, 
respectively,  and  are  far  too  large  for  the  algorithm  to  handle.  The  diffusion-based  boundary 
conditions  also  proved  more  difficult  to  obtain  a  stable  solution  for  in  the  model,  which  further 
contributed  to  the  mismatch  in  reconstructed  elasticity  contrast.  However,  the  results  obtained 
using  the  thin-plate  spline  method  are  encouraging  because  the  mean  error  was  0.43  voxel  units, 
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thereby  satisfying  the  threshold  (<0.5  voxel  units)  while  demonstrating  reconstruction  success. 
The  reconstruction  behavior  in  that  case  was  consistent  with  the  predicted  objective  function 
space  and  the  optimal  elasticity  contrast  was  found  to  be  within  6%  of  the  true  value.  This 
preliminary  result  appears  to  identify  the  use  of  thin-plate  spline  interpolation  as  a  strong 
candidate  for  generating  boundary  conditions  for  MIE.  The  implementation  for  that  method 
required  the  use  of  40  control  points,  which  is  seen  as  a  reasonable  choice  in  placing  fiducials  for 
data  acquisition  in  order  to  capture  the  extent  of  anticipated  deformation  processes.  A  more 
detailed  explanation  of  this  work  and  summary  of  experimental  results  is  provided  in  the 
accompanying  document  found  in  Appendix  E. 


Task  2(a)  stated:  “Design  and  construction  of  breast  phantoms  with  simulated  tumors.” 

We  have  opted  to  work  with  polyvinyl  alcohol  cryogel  (PVA-C)  because  of  its  relative  ease  and 
safety  in  manufacture  and  handling.  The  polymer  has  previously  been  described  as  being  a 
desirable  material  for  the  fabrication  of  both  ultrasound  and  MR  elastography  phantoms  for  use 
in  representing  the  brain,  blood  vessels,  and  breast  [11-13].  To  create  a  breast-shaped  fascimile, 
approximately  650  cc  of  8%  wt/vol  polymer  solution  is  molded  in  a  plastic  container  for  1  or  2 
freeze-thaw-cycles.  The  result  is  a  dome  shape  approximately  10-11  cm  at  the  base  and  tapering 
over  a  depth  of  about  5-6  cm. 

Because  of  its  use  of  intensity-based  image  similarity  metrics,  the  MIE  methodology 
requires  image  texture  in  the  variation  of  intensity  values  within  the  domain  to  perform  a 
registration  assessment  based  on  similarity.  While  the  natural  structures  of  the  breast  present 
little  concern  in  resulting  in  a  homogeneous  distribution,  introducing  usable  patterns  into  a 
phantom  poses  some  technical  challenges  in  order  to  not  have  intensity  correlate  too  strongly 
with  structural  differences  and  to  avoid  altering  the  elastic  properties  by  creating  a  heterogeneous 
composite  material.  A  method  that  has  been  found  to  be  potentially  satisfactory  for  MIE  involves 
doping  a  separate  quantity  of  polymer  solution  with  an  imaging  contrast  agent.  Before  the 
phantom  fully  polymerizes  through  freezing,  a  hypodermic  needle  is  used  to  create  multiple 
tracts  of  the  contrast-enhanced  slurry,  thereby  creating  regions  of  different  intensity  while 
maintaining  a  nearly  homogeneous  bulk  material.  Figure  14  below  shows  a  few  slices  of  a  PVA- 
C  phantom  constructed  in  this  manner  using  iohexol  suspension  (Omnipaque®,  GE  Healthcare, 
Chalfont  St.  Giles,  UK)  and  scanned  by  a  clinical  CT  unit. 


Figure  3.  Selected  views  of  a  PVA-C  breast  phantom  using  iodine  injections  and  imaged  by  CT. 
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Experimentation  has  also  led  to  a  protocol  that  should  allow  for  the  creation  and 
implantation  of  a  simulated  tumor  in  the  phantom.  Prior  to  manufacture  of  the  main  breast 
phantom,  a  spherical  inclusion  can  be  made  (out  of  a  significantly  stiffer  piece  of  PVA-C)  about 
2.5-cm  in  diameter  by  using  at  least  two  freeze-thaw-cycles  on  the  liquid  polymer.  This  inclusion 
can  then  be  suspended  within  the  main  breast  mold  to  let  the  bulk  material  solidify  around  it. 


Task  2(b)  stated:  “Prototype  construction  of  compression  device  including  design 
modifications  to  optimize  compatibility  with  current  imaging  modalities.”  and  Task  2(c) 
stated:  “Conduct  breast  phantom  experiments  to  investigate  efficacy  of  MIE  reconstruction 
while  varying  parameters  (e.g.  degree  of  compression)  during  imaging.” 

Two  sets  of  compression  devices  have  been  constructed  and  and  found  to  be  compatible  in 
magnetic  resonance  (MR)  and  X-ray  computed  tomography  (CT)  imaging  systems.  Both 
devices  are  primarily  composed  of  a  clear  acrylic  material  and  make  use  of  neoprene  air  bladders 
to  deliver  a  compression  of  up  to  5  cm.  The  first  device  is  a  rectangular  Plexiglas  frame  that  traps 
the  phantom  in  at  least  two  directions  with  a  sliding  wall  and  fixed  wall  which  houses  the  air 
bladder.  This  unit  has  been  used  on  PVA-C  phantoms  fabricated  as  described  above.  Figure  4 
below  illustrates  the  use  of  the  device  (in  this  case,  in  CT);  the  specifics  of  the  experiment  are 
described  in  Appendix  F. 

A  prototype  compression  chamber  that  is  more  clinically  oriented  has  been  designed  to  fit 
into  the  chassis  of  a  Philips  Intera  breast  coil  unit.  In  order  to  retain  the  geometry  of  the  overall 
device,  which  is  made  to  accommodate  95%  of  the  population  for  imaging  purposes  according  to 
the  manufacturer,  the  chamber  was  fabricated  from  by  cutting  cylindrical  segments  at  ~27°.  The 
air  bladders  are  attached  using  polycarbonate  pins  to  face  the  cranial  and  caudal  (supra-  and 
infra-mammary  aspects,  respectively)  surfaces  of  the  breast.  The  overall  assembly  is  then 
covered  with  an  expandable  nylon  sheet  (see  Figure  5). 

Image  acquisition  studies  in  conjunction  with  the  Vanderbilt  University  Institute  of 
Imaging  Science  have  recently  begun  to  demonstrate  the  clinical  feasibility  of  the  design. 

A  Philips  Achieva  3T  equipped  with  field  strengths  up  to  80  mT/m,  and  slew-rates  up  to 
200T/m/s  and  an  eight  channel  SENSE  breast  coil  has  been  used  for  all  breast  imaging.  The 
target  and  deformed  images  were  acquired  with  a  3D  Ti-weighted  high  resolution  isotropic 
volume  exam  (THRIVE)  that  includes  a  fat-nulling  inversion  pulse.  Pulse  sequence  parameters 
were  TR/TE=6.19ms/3.2ms  with  a  flip  angle  of  10°  and  a  NEX=1.  129  2-mm  thick  slices  were 
acquired  with  an  acquisition  matrix  was  400x400  zero  reconstructed  to  512x512  over  a  field  of 
view  of  20  cm2  to  25.6  cm2  (depending  on  breast  size).  This  sequence  design  is  currently  being 
utilized  for  initial  MIE  analysis  and  will  be  altered  accordingly  if  necessary.  Figure  6  below 
shows  a  sample  image  slice  of  a  patient  breast  with  a  tumor  obtained  from  the  MIE  clinical 
compression  chamber. 
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Figure  4.  Experimental  system  for  image  data  acquisition  on  phantoms.  A  polyvinyl  alcohol  cryogel  is  placed 
within  a  Plexiglas  chamber  with  its  surfaces  held  in  place  against  the  walls,  (a)-(c)  volumetric  segmentations  of  the 
CT  image  volume  taken  of  the  phantom  at  0%,  50%,  and  100%  inflation  of  the  air  bladder,  (d)  photograph  of  the 
setup.  Compression  is  delivered  through  an  air  bladder  (right  wall  in  picture)  inflated  manually  through  a  bulb 
adapted  from  a  standard  sphygmomanometer,  [from  Appendix  F] 
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Figure  5.  Left:  Photograph  of  assembly  looking  down  into  the  Philips  Intera  breast  coil  unit.  Right:  front  view  of 
the  device,  which  is  clamped  to  the  chassis  but  does  not  interfere  with  the  imaging  field  of  view. 


Figure  6.  Example  slices  of  an  MR  image  obtained  of  a  breast  (a)  before  and  (b)  after  compression  as  delivered  by 
the  clinical  MIE  chamber. 


Key  Research  Accomplishments 

Although  the  basic  requirements  of  a  static  compression  and  pre-/post-  image  analysis  that  are 
the  foundation  of  MIE  are  simple  conceptually  compared  with  other  elastography  methods,  the 
implementation  of  a  robust  system  is  a  significant  engineering  challenge.  The  flexibility  of  the 
parallelized  code  base  is  for  this  project  has  provided  an  interesting  avenue  of  analysis  of  a  large- 
scale  nonlinear  optimization  problem.  In  addition,  the  relatively  easy  translation  of  data 
acquisition  system  designs  to  the  clinic  is  a  key  breakthrough  for  the  project  to  demonstrate  the 
safety  and  applicability  of  the  method  in  a  real-world  setting,  whether  on  phantom  or  human 
subjects. 
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Reportable  Outcomes 


Work  on  the  MIE  method  has  so  far  resulted  in  multiple  conference  papers.  Prior  work  that  was 
completed  in  the  reporting  timeline  resulted  in  a  peer-reviewed  journal  publication.  These  items 
have  provided  the  foundation  for  a  thesis  proposal  approved  by  the  graduate  school,  as  well  as  a 
forthcoming  manuscript  currently  being  prepared. 

Poster  presentations 

Vanderbilt  University  Medical  Scientist  Training  Program  retreat  (July  2006).  Brentwood,  TN. 
Conference  papers 

Ou  JJ,  Barnes  SL,  and  Miga  MI,  "Application  of  multi-resolution  modality  independent 
elastography  for  detection  of  multiple  anomalous  objects,"  in  Medical  Imaging  2006: 

Physiology,  Function  and  Structure  from  Medical  Images,  San  Diego,  CA,  Feb  2006,  pp. 
614310-1  to  614310-9. 

Ou  JJ,  Barnes  SL,  Miga  MI,  "Preliminary  testing  of  sensitivity  to  input  data  quality  in  an 
elastographic  reconstruction  method,"  in  IEEE  International  Symposium  on  Biomedical  Imaging, 
Arlington,  VA,  Apr  2006,  pp.  948-951. 

Schuler  DR,  Ou  JJ,  Barnes  SL,  Miga  MI,  “Automatic  surface  correspondence  methods  for  a 
deformed  breast,”  in  Medical  Imaging  2006:  Visualization  and  Image-Guided  Procedures,  San 
Diego,  CA,  Feb  2006,  pp.  614125-1  to  614125-8. 

Ou  JJ,  Ong  RE,  Miga  MI,  “An  Evaluation  of  3D  Modality  Independent  Elastography 
Robustness  to  Boundary  Condition  Noise”,  SPIE  Medical  Imaging  2007.  Accepted  Oct  2006,  In 
press. 

Ong  RE,  Ou  JJ,  Miga  MI,  “Using  Laplace’s  equation  for  non-rigid  registration  of  breast 
surfaces”,  SPIE  Medical  Imaging  2007.  Accepted  Oct  2006,  In  press. 

Miga  MI,  Ou  JJ,  Ellis  DL,  “An  elastography  framework  for  use  in  dermoscopy”,  SPIE  Medical 
Imaging  2007.  Accepted  Oct  2006,  In  press. 


Conclusions 

The  current  results  and  progress  denoted  in  this  report  are  within  the  proposed  statement  of  work 
and  are  encouraging  towards  completion  of  the  overall  objectives  with  further  effort.  No 
significant  deviations  are  reported  at  this  time. 
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ABSTRACT 

This  work  extends  a  recently  realized  inverse  problem  technique  of  extracting  soft  tissue  elasticity  information  via  non- 
rigid  model-based  image  registration.  The  algorithm  uses  the  elastic  properties  of  the  tissue  in  a  biomechanical  model  to 
achieve  maximal  similarity  between  image  data  acquired  under  different  states  of  loading.  A  new  multi-resolution,  non¬ 
linear  optimization  framework  has  been  employed  which  allows  for  improved  performance  and  object  detection.  Prior 
studies  have  demonstrated  successful  reconstructions  from  images  of  a  tissue-like  thin  membrane  phantom  with  a  single 
embedded  inclusion  that  was  significantly  stiffer  than  its  surroundings.  For  this  investigation,  a  similar  phantom  was 
fabricated  with  two  stiff  inclusions  to  test  the  effectiveness  of  this  method  in  discriminating  multiple  smaller  objects. 
Elasticity  values  generated  from  both  simulation  and  real  data  testing  scenarios  provided  sufficient  contrast  for  detection 
and  good  quantitative  localization  of  the  inclusion  areas. 

Keywords:  Elastography,  elasticity  imaging,  multi-resolution  methods,  image  similarity,  finite  elements 

1.  INTRODUCTION 

The  practice  of  palpating  soft  tissue  structures  in  the  course  of  the  clinical  physical  exam  has  had  a  long-standing 
history  of  providing  correlation  of  improper  stiffness  with  pathology.  The  ability  to  characterize  the  mechanical 
properties  of  tissue  is  a  potential  source  of  additional  information  relevant  for  detection  and  diagnosis  of  a  disease 
process,  and  has  implications  for  the  assessment  of  treatment.  One  way  in  which  this  could  be  achieved  in  a  minimally 
invasive  manner  is  by  analyzing  tissue  deformation  through  imaging  and/or  image  processing  techniques,  which  is  a 
central  goal  of  the  field  of  elastography  [1].  Application  of  such  methods  to  the  interrogation  of  the  breast  [2,3],  skin 
[4-6],  prostate  [7],  and  other  accessible  organ  systems  is  an  emerging  area  of  research. 

Many  of  the  current  elastography  methods  are  founded  in  ultrasound  (US)  and  magnetic  resonance  imaging 
(MR)  and  involve  the  estimation  of  induced  displacements  within  the  tissue  of  interest  to  infer  the  elasticity  distribution. 
We  have  pursued  the  development  of  a  reconstruction  method  utilizing  quasi-static  deformation  and  image  similarity 
metrics  that  has  been  termed  ’modality-independent  elastography’  (MIE)  [8-10]  because  of  its  potential  to  handle  native 
anatomical  image  data  from  different  modalities  with  simple  modification  to  the  acquisition  procedure.  Common 
problems  facing  all  of  these  methods  involve  limitations  with  the  accurate  recovery  of  elastic  property  values,  detection 
of  small  lesions  in  tissue,  and  the  resolution  of  multiple  discrete  lesions  [11,12].  Building  upon  recent  study  involving  a 
single  focal  lesion  [6],  the  objectives  of  this  work  were  to  challenge  the  ability  of  the  MIE  method  to  reconstruct  a 
scenario  of  two  small  inclusions  embedded  in  a  homogeneous  domain  and  to  further  explore  the  feasibility  of  the 
method  in  handling  image  data  from  different  imaging  modalities.  This  was  accomplished  by  performing  simulated 
reconstructions  using  images  obtained  from  X-ray  computed  tomography  (CT),  MR,  and  digital  photography  and  then  a 
reconstruction  from  a  real-world  experiment  using  a  thin  phantom  membrane. 

2.  METHODS 


2.1  Elastographic  reconstruction  framework 

The  conceptual  framework  for  our  elastographic  reconstruction  has  been  previously  described  in  [6,8-10].  In  brief,  an 
image  of  a  tissue  of  interest  ( source )  is  deformed  by  a  biomechanical  computer  model  and  compared  against  an 
acquired  image  of  the  same  tissue  in  a  mechanically  loaded  state  {target).  The  deformation  and  comparison  is  repeated 
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using  systemic  updates  of  elasticity  parameters  until  a  suitable  match  in  intramodal  image  similarity  is  achieved  in  a 
least  squares  manner  to  satisfy  a  multi-resolution,  non-linear  optimization  scheme.  This  process  can  be  classified  as  an 
inverse  problem,  with  model-based  deformation  of  the  source  image  representing  the  forward  problem.  Each  of  the 
three  major  components  (model,  image  comparison,  and  optimization)  is  described  in  more  detail  in  the  following 
sections,  and  a  flow  chart  representation  of  the  overall  process  is  included  in  Figure  1. 

2.1.1  Biomechanical  model 

A  central  component  to  the  model-based  inverse  problem  is  the  manner  in  which  the  continuum  is  represented.  While 
the  constitutive  model  that  best  describes  tissue  deformation  mechanics  is  more  complex,  for  this  work,  linear  isotropic 
elasticity  has  been  employed.  The  partial  differential  equation  that  expresses  a  state  of  mechanical  equilibrium  can  be 
written  as  [13]: 

V  •  cr  =  0  (1) 


where  cris  the  Cartesian  stress  tensor. 

For  the  purposes  of  the  following  experimentation,  we  also  apply  either  the  plane  stress  or  plane  strain 
approximations  to  the  thin  membrane  and  breast  cross-section  trials,  respectively.  The  direct  consequence  of  this  is  a 
reduction  of  the  36  stiffness  constraints  in  the  general  3D  formulation  of  Cauchy’s  Faw  to  the  two  parameters  of 
Young’s  modulus  (E)  and  Poisson’s  ratio  (v)  in  2D.  These  simplifications,  while  significant,  are  appropriate 
descriptions  of  sufficiently  thin  and  thick  systems  under  planar  loading.  In  plane  stress, 
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describes  the  constitutive  relationship  between  the  Cartesian  stress  tensor  [gx,  ay,  xxy]  and  strain  tensor  [cx,  sy,  yxy]. 
Similarly,  in  plane  strain, 
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A  finite  element  (FE)  model  using  triangular  elements  is  constructed  from  the  source  image  and  assigned  appropriate 
boundary  conditions  based  on  estimated  displacement  or  stress  (i.e.  Dirichlet  and  Neumann  conditions,  respectively). 
The  standard  Galerkin  method  of  weighted  residuals  [14]  is  used  to  construct  and  solve  the  system. 

2.1.2  Image  deformation  and  comparison 

To  further  describe  the  reconstruction  process,  we  introduce  some  additional  terminology  at  this  point.  The  model 
domain  is  equivalent  to  the  total  area  of  the  FE  mesh  constructed  using  the  source  image  as  stated  above  and  contains 
the  relevant  elasticity  information.  The  model  domain  is  partitioned  by  a  K-means  clustering  of  the  element  centroids 
(MATFAB  R14,  Mathworks,  Natick,  MA)  into  N  number  of  regions ,  each  of  which  has  a  distinct  set  of  spatially 
homogeneous  elastic  properties.  Subdividing  in  this  manner  allows  for  the  implementation  of  the  multi-resolution 
reconstruction  whereby  progressively  finer  spatial  distributions  of  elasticity  parameters  are  utilized  in  the  process,  a 
method  that  improves  upon  previous  versions  using  only  a  single  resolution  [8-10].  Analogously,  the  comparison 
domain  is  an  area  specified  by  semi-automated  segmentation  on  the  target  image  and  contains  information  pertaining  to 
image  similarity.  The  comparison  domain  is  separated  into  M  number  of  rectangular  zones  containing  approximately 
equal  numbers  of  pixels. 
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The  reconstruction  algorithm  begins  by  assigning  an  initial  Young’s  modulus  value  to  each  of  the  regions  at 
the  coarsest  resolution.  Poisson’s  ratio  is  held  constant  at  v  =  0.485  to  represent  a  nearly  incompressible  material.  The 
FE  model  is  solved  to  determine  the  nodal  mesh  displacements,  which  are  in  turn  used  to  deform  the  source  image.  This 
model-deformed  image  is  then  compared  to  the  target  image  for  every  zone  using  an  intensity-based  image  similarity 
metric.  While  a  number  of  methods  are  available  for  such  a  task,  here,  we  utilize  the  correlation  coefficient  (CC)  [15] 
throughout,  as  it  has  empirically  demonstrated  superior  performance  over  other  metrics  such  as  the  sum  of  squared 
differences  and  normalized  mutual  information. 


2.1.3  Optimization  scheme 


Let  T  be  a  function  that  represents  the  model-based  image  deformation  and  takes  as  its  input  a  vector  of  elastic  modulus 
values  E  of  length  N  that  corresponds  to  the  current  distribution  of  regions  in  the  model  domain.  Then  for  two 
distributions  of  modulus  values  Ex  and  E2,  the  similarity  between  the  images  produced  by  J(Ei)  and  T( E2)  is  the  vector 
S  of  length  M  containing  evaluations  of  the  correlation  coefficient  corresponding  to  the  distribution  of  zones  in  the 
comparison  domain.  The  elasticity  parameter  optimization  can  be  written  as  the  minimization  of  the  least  squares  error 
objective  function 


vp  _  c  _  e 

1  y^TRUE  ° EST\ 
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where  STRUE  is  the  set  of  similarity  values  achieved  when  comparing  the  target  image  to  itself,  SEst  is  the  similarity 
between  the  model-deformed  source  and  the  target  images  using  current  estimates  of  the  elastic  modulus  distribution, 
and  |*|  denotes  the  vector  L2  norm.  By  definition,  STrue  is  the  maximum  value  for  the  similarity  metric  (max  CC  =1). 
Using  a  Levenberg-Marquardt  approach,  the  residual  form  of  equation  (4)  becomes 

[lTJ  +  all&E}  =  [jT -  SEST }  (5) 

where  J  =  dSESx/dE  is  the  Jacobian  matrix  of  size  MxN  and  I  is  the  N x  N  identity  matrix.  Because  JTJ  is  typically  an 
ill-conditioned  term,  the  regularization  parameter  a  is  determined  using  the  methods  described  in  [16].  Modulus  values 
of  the  regions  at  a  given  resolution  are  updated  by  AE  until  an  error  tolerance  is  reached  or  a  maximum  number  of 
iterations  have  been  completed.  Upon  reaching  a  stopping  criterion,  the  material  property  description  is  interpolated 
onto  the  next  (i.e.  finer)  resolution  and  the  above  steps  are  repeated.  Spatial  averaging  of  modulus  values  within  the 
model  domain  and  solution  relaxation  between  successive  resolution  levels  are  also  utilized  to  improve  the  stability  of 
the  optimization. 
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Figure  1.  Flow  chart  of  elastographic  reconstruction  framework. 

2.2  Reconstruction  experiments 

A  two-material  phantom  membrane  of  simulated  skin  had  been  previously  constructed  [6]  using  Smooth-On™ 
polyurethanes  (Smooth-On,  Easton,  PA)  designated  by  the  manufacturer  as  Evergreen  10  and  Evergreen  50.  These 
materials  have  essentially  indistinguishable  colors  but  vary  significantly  in  their  elastic  modulus  values,  so  the  former 
was  used  as  the  bulk  material  and  the  latter  for  stiff  objects.  From  material  testing,  the  elastic  modulus  contrast  was 
expected  to  be  approximately  5.7:1.  The  phantom  was  made  to  contain  two  circular  stiff  inclusions  1.5  cm  in  diameter 
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embedded  near  opposing  comers  of  a  rectangular  field  of  bulk  material  measuring  15  cm  x  14  cm.  A  black  permanent 
marker  was  used  to  place  a  pattern  of  regularly  spaced  (-1  cm)  grid  lines  across  the  membrane.  The  thin  membrane  was 
securely  clamped  along  two  opposite  edges  and  then  subjected  to  a  uniaxial  tensile  displacement  (-8%  strain)  by  means 
of  a  milling  vise.  A  commercial  webcam  (Logitech  QuickCam  Pro  4000,  960  x  1280  pixel  resolution)  was  rigidly 
mounted  above  the  membrane  to  acquire  image  pairs  of  the  pre-  and  post-stretched  states. 

To  initially  test  the  method  regarding  the  two-inclusion  scenario,  a  simulation  using  the  source  image  of  the 
membrane  was  performed  by  deforming  it  with  a  prescribed  model  (plane  stress)  of  known  boundary  displacements  and 
elasticity  parameters  to  generate  a  target  image;  high  modulus  values  were  assigned  to  elements  bounded  by  a 
segmentation  of  the  inclusion  locations.  A  reconstruction  was  then  performed  using  the  actual  image  data  acquired  as 
described  above.  In  both  cases,  resolutions  of  N=  16,  64,  256,  512,  and  800  regions  and  M=  400  zones  were  used.  The 
results  of  the  idealized  and  real  data  reconstructions  are  shown  in  Figures  4  and  5,  with  further  quantitative  evaluation  in 
Table  1. 

In  order  to  examine  the  robustness  of  the  method  regarding  its  use  of  data  from  differing  sources,  simulation 
reconstructions  were  performed  using  image  slices  extracted  from  breast  image  volumes  obtained  from  CT  and  MR 
scans  (see  Figure  3).  Although  these  were  taken  from  two  different  patients,  the  images  were  selected  to  be 
approximately  corresponding  slices  ~2  cm  away  from  the  chest  wall  in  the  coronal  orientation  of  the  standard 
anatomical  position.  The  simulations  were  set  up  in  the  same  manner  as  for  the  digital  photographs,  using  either  one  or 
two  inclusions  of  about  1  cm  in  diameter  embedded  within  the  true  elasticity  distribution  and  a  small  compression  (-8% 
strain)  in  the  cranial-caudal  direction.  The  relative  stiffness  of  the  inclusions  was  designated  to  be  5.7:1  for  consistency 
with  the  material  testing  data  and  also  because  the  value  is  fairly  representative  of  breast  tumor  properties  [17].  The 
plane  strain  model  approximation  was  used  in  the  breast  simulation  trials,  progressing  through  resolutions  of  N=  24,  64, 
256,  and  576  regions  using  M  =  200  zones.  The  reconstruction  method  was  then  run  for  all  four  test  cases,  and  the 
results  are  presented  in  Figures  6  and  7  and  Table  2. 


Figure  2.  (Left  to  right):  Phantom  membrane  in  undeformed  state  (source  image),  under  deformation  (target  image),  and  difference 
image.  Arrows  in  the  left  panel  indicate  the  positions  of  the  two  stiff  inclusions. 


2.3  Reconstruction  evaluation 

The  fidelity  of  the  elasticity  reconstruction  was  evaluated  on  its  ability  to  detect  the  presence  of  an  inclusion  based  on 
classification  of  the  material  property  distribution,  and  the  retrospective  accuracy  of  localizing  the  lesions.  The  elastic 
properties  as  a  whole  were  treated  as  a  Gaussian  mixture  of  two  classes  and  separated  by  a  threshold  established  via  the 
method  described  in  [18].  The  likelihood  of  detecting  a  lesion  in  the  elasticity  image  was  found  using  the  contrast- to- 
noise  ratio  as  defined  by  [12,19]: 


/2 (Ml  ~Maf 

V 


CNR  = 


(6) 


Figure  3.  Images  slices  of  breast  tissue  extracted  from  a  CT  volume  (left)  and  MR  volume  (right)  used  in  simulation  study  of  the 
ability  of  the  reconstruction  method  to  utilize  disparate  image  data  types. 


where  ju  and  o2  are  the  sample  mean  and  variance  of  a  material  property  distribution  and  the  subscripts  L  and  B  denote 
the  lesion  and  bulk  material  classes,  respectively.  As  a  quantitative  assessment  of  the  localization  of  the  lesion(s),  the 
positive  predictive  value  of  correctly  identifying  a  lesion  material  within  the  known  segmented  region  of  the  inclusions 
was  used  as  a  'quality  of  reconstruction  score'  (QRS).  This  value  is  significant  because  identification  of  the  lesion 
border  and  material  classification  are  done  independently,  so  any  user  knowledge  of  the  test  scenario  does  not  influence 
the  performance  of  the  measure.  Cutoffs  for  successful  detection  and  localization  were  set  at  CNR>2.2  as  noted  by  [12] 
and  QRS>80%  as  determined  by  prior  study  in  our  laboratory.  The  average  modulus  contrast  is  found  from  the  ratio  of 
the  means  of  the  two  material  classes,  and  a  peak  modulus  contrast  value  is  also  reported  by  taking  the  ratio  of  two 
manually  selected  homogeneous  regions  of  approximately  equal  area  known  to  be  representative  of  the  two  materials. 
It  should  be  noted  that  in  other  work  not  presented  here,  the  definition  of  QRS  included  a  weighting  factor  provided  by 
the  estimated  reconstruction  modulus  contrast,  but  for  the  current  purposes,  only  localization  accuracy  was  considered 
to  maintain  an  objective  evaluation  of  inclusion  detection. 


3.  RESULTS 

Figure  4  demonstrates  the  ability  of  the  reconstruction  method  to  produce  an  elasticity  map  from  the  simulation  data 
with  good  localization  of  the  inclusions  that  are  easily  visually  distinguishable  from  the  surrounding  bulk  material.  The 
progression  through  resolutions  of  N  =  64,  256,  512,  and  800  regions  shows  improving  delineation  of  the  inclusions  and 
elastic  contrast.  Figure  5  demonstrates  a  similar  behavior  for  the  reconstruction  of  the  acquired  phantom  membrane 
data,  with  both  spatial  definition  and  modulus  contrast  increasing  with  the  finer  discretization.  Table  1  summarizes  the 
quantitative  evaluation  of  the  reconstructions  in  both  simulation  and  phantom  trials,  including  CNR,  contrast  ratio,  and 
QRS  values.  The  CNR  values  are  sufficient  to  allow  for  discrimination  of  the  two  materials  and  the  identification  of  the 
inclusions  was  determined  to  be  accurate  in  both  cases.  The  reconstruction  of  the  phantom  membrane  does  show  some 
misclassification  along  the  border  where  the  deformation  was  applied  as  well  as  in  the  comer  adjacent  to  one  of  the 
inclusions  (see  Figure  5d). 

Figures  6  and  7  show  the  final  reconstruction  results  for  the  CT  and  MR  breast  slice  simulations  using  either 
one  or  two  inclusions.  In  both  test  scenarios,  the  resolvability  of  the  stiffer  material  was  found  to  be  adequate  according 
to  the  CNR  threshold,  but  definitely  higher  in  the  MR-derived  elasticity  images.  Localization  of  the  inclusions  yielded 
excellent  QRS  values  in  reconstmctions  using  either  modality,  again  higher  (though  slightly)  for  the  MR  images. 
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(a)  (b) 


(c)  (d) 


Contrast  Ratio 


Figure  4.  Reconstruction  of  the  simulated  membrane  deformation  using  idealized  model  parameters,  progressing  through  finer 
resolution  distributions  (a)-(d)  of  64,  256,  512,  and  800  regions. 


Figure  5.  Reconstruction  of  the  actual  membrane  data.  A  faint  contour  in  (d)  is  present  to  demarcate  the  actual  position  of  the  stiff 
inclusions.  Again,  panels  (a)-(d)  demonstrate  the  effect  of  the  multi-resolution  method  in  utilizing  64,  256,  512,  and  800  regions  to 
better  capture  the  shape  and  location  of  the  inclusions. 


Table  1.  Quantitative  reconstruction  evaluations. 


Avg  CR 

Max  CR 

CNR 

QRS  (%) 

Simulation 

2.7 

4.0 

4.4 

97.7 

Phantom 

2.6 

4.1 

2.8 

88.5 
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Figure  6.  Reconstructions  of  simulation  trials  for  the  CT  breast  slice  using  a  single  inclusion  (left)  and  two  inclusions  (right).  The 
true  inclusion  boundaries  are  overlaid  in  each  elasticity  image. 


Figure  7.  Reconstructions  (bottom  row)  of  simulation  trials  for  the  MR  breast  slice  using  a  single  inclusion  (left)  and  two  inclusions 
(right).  The  true  elasticity  distributions  are  also  shown  (top  row)  for  comparison. 


Table  2.  Quantitative  reconstruction  evaluations. 


Avg  CR 

Max  CR 

CNR 

QRS  (%) 

CT  (1  inclusion) 

2.1 

3.1 

3.0 

97.6 

CT  (2  inclusions) 

2.0 

2.6 

3.5 

96.9 

MR  (1  inclusion) 

2.8 

3.7 

20.0 

100 

MR  (2  inclusions) 

2.7 

3.7 

5.7 

99.8 
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4.  DISCUSSION 

The  results  of  the  phantom  membrane  experiment  are  encouraging  because  of  their  similarity  to  the  idealized 
simulation.  Despite  nonlinear  model-data  mismatch,  out-of-plane  distortions  during  stretching,  and  possible  boundary 
condition  inaccuracies,  the  elasticity  reconstruction  demonstrated  good  localization  of  the  two  small  inclusions.  The 
majority  of  the  problems  in  reconstruction  are  mostly  likely  due  to  noise  incurred  in  the  mapping  of  the  boundary 
displacements.  It  should  be  noted  that  the  phantom  reconstruction  was  achieved  with  a  non-pigmented  lesion  (see 
Figure  2,  arrows),  indicating  that  deflections  of  the  image  structure  are  capable  of  driving  the  image  similarity  metric  of 
the  reconstruction  process.  This  does  intuitively  suggest  that  some  metric  for  rating  the  complexity  and  density  of 
image  pattern  in  relation  to  algorithm  success  may  be  important  and  is  currently  under  investigation.  Preliminary  data 
not  presented  in  this  work  indicates  that  such  a  threshold  does  exist  for  image  data  that  can  be  properly  analyzed  by  the 
current  framework.  The  modality  independence  of  the  method  is  also  supported  by  the  results  here;  clearly,  the 
Hounsfield  units  of  CT,  floating  point  values  from  an  MR  volume,  and  the  luminance  captured  by  the  CCD  sensor  of  a 
digital  camera  are  quite  different  types  of  data  to  handle  because  they  are  based  on  different  physical  principles.  The 
simulation  reconstructions  demonstrate  that  the  method  is  indifferent  to  these  differences  by  treating  the  data  as  an 
arbitrary  range  of  intensities  and  will  converge  towards  the  true  elasticity  distribution  based  on  the  image  pattern 
available.  This  is  a  possible  explanation  for  the  qualitatively  more  satisfactory  results  from  the  MR  simulations 
compared  to  the  CT  trials  because  the  distribution  of  intensities  from  the  former  modality  yielded  a  more  diversified 
histogram,  an  attribute  that  should  naturally  aid  an  intensity-based  metric. 

While  an  ideal  reconstruction  would  also  be  accurate  in  characterizing  a  lesion  by  its  modulus  contrast,  our 
focus  in  the  study  was  to  test  the  ability  of  the  method  to  detect  and  localize  the  inclusions.  In  previous  experimentation 
with  reconstructions  of  single  focal  lesions,  we  have  been  generally  successful  in  achieving  a  contrast  ratio  within  25% 
of  the  true/expected  value.  It  is  somewhat  troubling  that  the  contrast  ratios  calculated  here  did  not  meet  that  criterion, 
although  the  experiments  with  the  phantom  membrane  came  fairly  close  (28%).  However,  these  results  underscore  the 
difficulty  of  the  scenarios  in  not  only  having  to  deal  with  multiple  inclusions  but  quite  small  ones  in  both  the  true 
physical  sense  and  also  the  scale  of  the  domain.  Any  of  the  given  inclusions  tested  in  simulation  and  with  the  real  data 
were  detected  within  a  homogeneous  domain  approximately  an  order  of  magnitude  larger  (e.g.,  1.5-cm  lesions  in  a  15 
cm  x  14  cm  domain  for  the  phantom).  The  expectation  of  being  able  to  identify  with  any  confidence  the  presence  of  the 
inclusion  is  comparable  to  the  observations  made  in  [12]  where  the  test  of  finding  a  single  5 -mm  lesion  within  a  4  cm  x 
5  cm  domain  proved  to  be  the  most  problematic.  Therefore,  the  localization  of  the  lesions  as  determined  by  the  CNR 
and  QRS  metrics  is  deemed  to  be  a  success,  and  further  investigation  into  the  nature  of  the  method  with  respect  to  the 
scale  of  the  lesion  and  domain  is  warranted. 


5.  CONCLUSIONS 

In  this  work,  we  have  presented  further  testing  of  a  method  for  recovering  elasticity  parameters  by  maximizing  the 
similarity  between  images  of  a  tissue  of  interest  acquired  under  two  different  states  of  quasi-static  loading  within  the 
context  of  an  inverse  problem.  The  specific  experiments  presented  here  examined  the  effectiveness  of  the  technique  for 
the  detection  of  multiple  small  discrete  focal  lesions  embedded  in  an  otherwise  homogeneous  medium,  as  well  as 
further  proof-of-concept  work  in  its  applicability  to  utilize  image  data  from  various  modalities.  In  both  cases,  the 
method  provided  accurate  localization  of  the  lesions  based  on  the  reconstruction  of  relevant  elasticity  contrast.  Because 
the  biomechanical  model,  multi-resolution  optimization,  and  image  acquisition  are  each  modular  components  of  the 
framework,  this  elastographic  reconstruction  technique  is  readily  extensible  for  added  sophistication,  and  there  is 
ongoing  work  to  enhance  the  methodology  with  more  complex  models  and  advances  in  imaging  technology. 
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Modality  Independent  Elastography 
(MIE)  Concepts 


•  (Solid)  tumors  are  usually  stiffer  than  surrounding 
tissue 


•  Soft  tissue  interrogation  of  various  organ  systems 
(e.g.  skin,  liver,  prostate,  breast)  for  tumor 
detection 


•  Elastography  gives  representation  of  a  structure 
according  to  its  mechanical  properties 

•  Deformation  processes  indicative  of  material 
inhomogeneity  can  be  captured  by  imaging  and 
approximated  with  modeling 

•  Associate  form  and  function  through  image  analysis 
separate  from  modality  acquisition 
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MIE  Components 

•  (1)  Biomechanical  FE  model  of 

soft-tissue  deformation 

•  Conservation  of  stresses  (continuum) 

•  Constitutive  stress-strain  relation  (Hooke's  Law) 


<j  =  Es 


MIE  Components  (cont.) 

•  (2)  Similarity  measure  for 

comparing  images 

•  Acquired  "pre-"  (source)  &"post-"  (target) 
quasi-static  deformation 

•  Intensity-based  registration  metrics 
•  MI,  NMI,  SSD,  CC,  GC 
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MIE  Components  (cont.) 


•  (3a)  Optimization  routine  to  update 
material  properties  in  the  model 


•  Objective  function  based  on  similarity 


MIE  Components  (cont.) 

•  (3b)  Discretization  of  elasticity 
distribution  and  image  data 

•  Multi-resolution  K-means  clustering  of  elements 
("regions") 

•  Sampling  of  image  comparison  area  ("zones") 
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MIE  Framework 


Medical  Physics,  vol.  32,  no.  5,  pp.  1308-1320,  2005 
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Study  Objectives: 

Further  Testing  of  MIE 

•  Modality  independence 

Digital  photography 

X-ray  computed  tomography  (CT) 

Magnetic  resonance  (MR) 

•  Two  (small)  inclusions 

•  Simulation  and  phantom  membrane 
study 
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Evaluating  MIE 

•  Classify  reconstruction 

Two-class  Gaussian  mixture  model 

•  Detectability  via  elasticity  image  contrast 


•  Localization  accuracy 

Positive  predictive  value  of  identifying  lesion 
material  in  correct  location 


CT  breast  slice  -  simulation 

Source  Targetl  Target2 
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MR  breast  slice  -  simulation 

Source  Targetl  Target2 
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Summary 

•  Modality  independence  via  simulation  for 
handling  various  data  types 

•  Multi-resolution  approach  potentially 
improves  optimization  convergence 

•  Two  small  stiff  inclusions  reconstructed  in 
phantom  membrane  experiment 

•  Detectability  accomplished  via  CNR 

•  Localization  successful  as  evaluated  by  QRS 
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<finis> 


Future  Directions  /  Research  Questions 

•  Biological  tissues  are  not  typically  linearly 
elastic 

•  Need  for  accurate  boundary  conditions 
creates  dependence  on  segmentation 
methodology 

•  Not  all  data  sets  necessarily  contain 
sufficient  information  for  elastographic 
reconstruction 
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ABSTRACT 

An  elastographic  reconstruction  method  has  been  developed 
to  recover  the  material  properties  of  soft  tissue  by  model- 
based  analysis  of  image  data  acquired  at  different  states  of 
mechanical  loading.  The  algorithm  utilizes  image  similarity 
as  part  of  the  cost  function  for  a  multi-resolution,  non-linear 
optimization.  Previous  work  with  a  phantom  membrane 
used  for  simulated  dermoscopic  application  has  prompted 
this  preliminary  investigation  of  the  relative  effects  of 
additive  image  noise  and  boundary  condition  determination 
errors  on  the  performance  of  the  method.  The  results  as 
quantified  by  elasticity  contrast  and  localization  accuracy 
indicate  that  the  reconstruction  process  is  robust  in  the 
presence  of  realistic  levels  of  image  corruption  and  tolerates 
the  majority  of  boundary  condition  mapping  errors. 


1.  INTRODUCTION 

The  practice  of  palpating  soft  tissue  structures  in  the  course 
of  the  physical  exam  for  assessing  tissue  health  has  had  a 
long-standing  clinical  history  of  providing  correlation 
between  improper  stiffness  and  pathology.  The  ability  to 
characterize  the  mechanical  properties  of  tissue  is  therefore 
a  potential  source  of  information  relevant  for  both  diagnosis 
and  prognosis.  One  way  in  which  this  could  be  achieved  in  a 
non-invasive  manner  is  through  analysis  of  tissue 
deformation  with  imaging  and  image  processing  techniques, 
which  is  a  central  goal  of  the  field  of  elastography  [1]. 

The  conceptual  framework  for  our  elastographic 
reconstruction  has  been  previously  described  in  [2-4].  In 
brief,  images  of  a  tissue  of  interest  are  acquired  in  an  initial 
{source)  and  then  mechanically  loaded  state  {target).  The 
source  image  is  deformed  by  a  prescribed  computational 
model  and  compared  to  the  target.  This  is  repeated  in  an 
iterative  process  using  updates  to  the  elasticity  parameters  of 
the  model  as  generated  by  a  multi-resolution,  non-linear 
optimization  scheme  in  order  to  achieve  a  suitable  match  in 
image  similarity.  Because  the  goal  of  the  reconstruction  is  to 


determine  a  spatial  mapping  of  tissue  elasticity,  this  process 
can  also  be  classified  as  an  inverse  problem. 

Our  observations  during  the  ongoing  development  and 
testing  of  this  method  have  prompted  questions  concerning 
the  quality  of  data  necessary  and  sufficient  to  achieve 
satisfactory  results  (i.e.  fidelity  of  the  reconstructed 
elasticity  image).  The  primary  inputs  to  the  reconstruction 
method  are  the  acquired  images  and  the  delineated  boundary 
conditions  on  the  region  of  interest.  While  it  is  clearly 
preferable  to  have  idealized  data,  in  reality,  both  inputs 
involve  varying  levels  of  manual  interaction.  As  an  initial 
study,  we  have  sought  to  test  the  effects  of  degradation  in 
data  quality  on  the  end  reconstruction  by  using  additive 
image  noise  and  randomized  boundary  condition  selection 
error. 

2.  METHODS 

2.1.  Elastographic  Reconstruction  Framework 

There  are  three  major  components  in  the  reconstruction 
framework:  a  biomechanical  model  of  tissue  response  to 
applied  deformation,  a  method  of  image  comparison,  and  an 
optimization  scheme.  For  the  current  version,  a  continuum- 
based  model  of  mechanical  equilibrium  using  isotropic 
Hookean  linear  elasticity  with  a  plane  stress  approximation 
is  employed  [5].  This  allows  for  a  reduction  of  the  general 
3D  formulation  of  Cauchy’s  Law  to  the  two  parameters  of 
Young’s  modulus  and  Poisson’s  ratio  in  2D.  The 
displacement  solution  of  the  finite  element  representation  of 
the  model,  solved  using  the  standard  Galerkin  method  of 
weighted  residuals  [6],  is  then  applied  to  the  nodes  of  a 
simple  triangular  mesh  based  on  the  source  image  domain  in 
order  to  perform  image  deformation.  The  mesh  is 
partitioned  by  K-means  clustering  (MATLAB  R14, 
Mathworks,  Nattuck,  MA)  into  N  number  of  regions ,  each 
of  which  describes  a  distinct  set  of  homogeneous  elastic 
properties  for  a  grouping  of  adjacent  elements.  This  allows 
for  implementation  of  the  multi-resolution  approach  by 
creating  a  hierarchy  of  increasingly  finer  spatial 
distributions  of  elasticity  parameters,  which  has  been  shown 
to  be  an  improvement  upon  previous  versions  using  only  a 
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single  resolution  [2,3].  A  second  discretization  is  performed 
to  divide  the  target  image  into  M  number  of  rectangular 
zones  containing  approximately  equal  numbers  of  pixels. 
The  deformed  source  image  is  compared  to  the  target  using 
an  intensity-based  image  similarity  metric  (here,  the 
correlation  coefficient  [7])  in  the  evaluation  of  the  least 
squares  error  objective  function 

M 

(S  -s  Y  t1) 

\^TRUE  °EST ) 

m= 1 

where  STrue  is  an  Mxl  vector  of  the  (maximum)  similarity 
values  achieved  when  comparing  the  target  image  to  itself 
and  Sest  is  the  Mxl  vector  of  similarity  between  the  target 
and  model-deformed  source  image  created  using  current 
estimates  of  the  elastic  modulus  distribution.  It  should  be 
noted  that  Sjrue  has  by  definition  a  value  of  1  for  the 
correlation  coefficient. 

The  minimization  of  equation  (1)  using  a  Levenberg- 
Marquardt  approach  takes  the  form 

[jTJ  +  cd]lAE}=  [/7  -  SEST }  (2) 

where  J  is  the  Jacobian  matrix  of  size  MxN  estimating 
dS/dE,  A E  is  the  Nxl  vectors  of  updates  to  the  current 
elasticity  values,  and  a  is  the  scalar  regularization  term  for 
the  Hessian  matrix  as  described  in  [8]. 

2.2.  Material  Preparation  and  Image  Acquisition 

For  our  simulation  purposes,  a  two-material  skin  phantom 
had  been  previously  constructed  [2]  as  a  thin  membrane 
measuring  15  cm  x  15  cm,  with  a  single  5 -cm  circular  stiff 
inclusion  embedded  in  the  center  (Figure  1).  The  phantom 
was  manufactured  with  Smooth-On™  polyurethanes 
(Smooth-On,  Easton,  PA)  Evergreen  10  and  Evergreen  50. 
These  materials  have  essentially  indistinguishable  colors  but 
vary  significantly  in  their  elastic  modulus  values,  so  the 
former  was  used  as  the  bulk  material  and  the  latter  for  the 
inclusion.  Based  on  material  testing,  the  expected  contrast 
ratio  of  Young’s  modulus  values  was  determined  to  be 
approximately  5.7:1.  A  black  permanent  marker  was  used 
to  place  a  pattern  of  regularly  spaced  (-1  cm)  grid  lines  on 
the  membrane.  The  membrane  was  clamped  along  two 
opposite  edges  and  then  stretched  in  a  uniaxial  fashion  by 
means  of  a  milling  vise.  A  commercial  webcam  (Logitech 
QuickCam  Pro  4000)  was  mounted  above  the  assembly  to 
acquire  image  pairs  of  the  membrane  in  pre-  and  post- 
stretched  states  (960  x  1280  pixel  resolution,  8-bit 
grayscale). 

2.3.  Reconstruction  Experiments 

Based  on  prior  work,  a  data  set  consisting  of  a 
particular  image  pair  and  associated  boundary  conditions 


Figure  1.  Experimental  phantom  membrane  system  (left)  and 
input  image  with  overlaid  finite  element  mesh  (right).  The 
inclusion  location  is  indicated  by  the  arrow  and  dotted  line.  The 
mesh  designates  the  actual  region  reconstructed. 

known  to  produce  a  satisfactory  reconstruction  was 
designated  as  the  gold  standard  for  the  remainder  of  the 
experiments  (Figure  1).  In  order  to  test  the  effect  of 
increasing  amounts  of  additive  noise  on  the  reconstruction 
algorithm,  Gaussian  random  fields  of  1,  5,  10,  15,  20,  25, 
and  30%  noise  were  applied  to  the  base  target  image  in  three 
separate  trials.  This  presents  a  challenge  that  ascertains  the 
ability  of  the  similarity  metric  and  objective  function  to 
discern  a  proper  match. 

The  current  method  for  selecting  Dirichlet  boundary 
conditions  on  the  finite  element  mesh  is  semi-automated  and 
requires  the  user  to  make  a  final  determination  on  point 
correspondence.  The  second  experiment  was  intended  to 
simulate  the  targeting  error  of  the  user  (e.g.  visual  cues  and 
input  device  control).  Each  test  involved  applying 
randomized  vectors  of  equal  magnitude  to  alter  the 
boundary  conditions  of  the  gold  standard  data  set.  Errors  of 
0.1,  0.2,  0.3,  0.5,  0.75,  1.0,  1.5,  and  2.0  mesh  units  (scaled 
to  be  equivalent  to  pixel  coordinates)  were  used  in  two 
separate  trials  for  a  total  of  16  reconstructions.  Sub-pixel 
magnitudes  were  included  after  determining  that  the 
accuracy  of  selecting  a  feature  point  in  the  image/mesh  was 
typically  less  than  or  equal  to  0.5  units  for  users  ranging 
from  moderate  to  expert  skill. 

For  all  reconstructions,  resolutions  progressing  through 
N  =  16,  36,  64,  144,  256,  and  400  regions  and  M  =  9 
similarity  zones  were  used;  domains  were  initialized  to 
homogeneous  elasticity  and  Poisson’s  ratio  held  constant  at 
0.485  to  represent  nearly  incompressible  mater ial(s). 

2.4.  Reconstruction  Analysis 

The  final  reconstructed  elasticity  values  were  modeled  as  a 
mixture  of  two  Gaussian  distributions,  and  a  threshold  was 
established  to  maximize  inter-class  variation  [9]  and 
subsequently  classify  each  region  as  bulk  or  stiff  material. 
Because  Dirichlet  boundary  conditions  are  exclusively  used 
in  these  reconstructions,  the  method  is  only  sensitive  to 
relative  differences  in  elasticity.  The  quantities  used  in 
evaluating  reconstruction  success  are  the  elasticity  contrast 
ratio,  localization  accuracy  of  the  inclusion,  and  an  overall 
measure  designated  the  ‘quality  of  reconstruction  score’ 
(QRS).  The  elasticity  contrast  ratio  (CR)  was  calculated 
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Figure  2.  Representative  reconstructions  with  image  noise.  From 
top  left:  1,  5,  10,  20,  25,  and  30%  additive  Gaussian  noise.  The 
reconstructions  are  visualized  as  two  materials,  with  whiter  areas 
indicating  higher  elasticity  contrast  values. 


Table  1.  Reconstruction  quality  under  noise  conditions 


Additive  image  noise 

%  Noise 

1 

5 

10 

15 

20 

25 

30 

LA 

0.92 

0.90 

0.91 

0.70 

0.69 

0.66 

0.56 

CR 

3.56 

3.45 

3.45 

3.24 

2.88 

2.83 

2.68 

Gold  standard:  LA  -  0.95,  CR  -  3.60 


Boundary  condition  error 

Err 

0.1 

0.2 

0.3 

0.5 

0.75 

1.0 

1.5 

2.0 

AE 

0.96 

3.32 

2.21 

102 

0.93 

32.2 

12.6 

7.66 

LA 

0.87 

0.92 

0.88 

0.59 

0.94 

0.86 

0.86 

0.96 

CR 

3.63 

3.68 

3.44 

2.91 

3.46 

3.71 

3.78 

3.30 

CR  = 
AE  = 

elasticity  contrast  ratio,  LA  =  localization  accuracy 
initial  alignment  error  (%),  Err  =  error  magnitude. 

from  the  mean  values  of  the  two  material  classes,  and  the 
positive  predictive  value  of  identifying  stiff  material  within 
the  independently  segmented  boundary  of  the  inclusion 
gives  a  measure  of  localization  accuracy  (LA).  The  quality 
of  reconstruction  is  simply  then  the  product  QRS  =  CR*LA, 
which  allows  the  user  to  consider  the  other  two  measures  in 
conjunction. 


3.  RESULTS 

Figures  2  and  3  show  examples  of  reconstructions  achieved 
under  various  image  noise  and  boundary  condition  errors, 
and  individual  localization  errors  and  contrast  ratio  values 
are  listed  in  Table  1.  Note  that  the  data  for  the  image  noise 
experiment  was  averaged  from  the  three  trials,  and  that  the 
data  presented  for  the  boundary  condition  experiment  is 
from  one  [representative]  trial.  Figure  4  is  a  plot  of  the 
reconstruction  quality  decreasing  with  increasing  image 
noise,  and  Figure  5  shows  the  reconstruction  quality  trend 
plotted  against  the  change  in  initial  alignment  error  (detailed 
in  the  following  section)  relative  to  that  of  the  gold  standard. 


□  DD 

□  DD 


Figure  3.  Representative  reconstructions  with  boundary  condition 
error.  Left  to  right:  0.1,  0.2,  0.3  units  (top  row);  0.75,  1.0,  2.0  units 
(middle  row,  trial  #1);  0.75,  1.0,  2.0  units  (bottom  row,  trial  #2). 
Error  magnitudes  greater  than  or  equal  to  0.5  mesh  units  are  not 
accurate  predictors  of  reconstruction  quality. 


4.  DISCUSSION 

From  visual  inspection  of  Figure  2,  it  is  apparent  that  the 
achieved  reconstruction  becomes  more  inaccurate  with 
increased  image  noise.  However,  the  ability  to  identify  and 
localize  the  stiff  inclusion  is  not  significantly  impaired  until 
a  noise  field  of  greater  than  10%  is  applied.  The  threshold 
was  found  by  determining  which  level  of  noise  provided  the 
best  minimum  sum  squared  error  fit  of  two  lines  to  the 
distribution  of  reconstruction  quality  in  Figure  4.  This 
would  indicate  that  the  similarity  metric  and  objective 
function  are  robust  to  intensity  deviations  of  about  6  gray 
levels.  While  Gaussian  noise  is  one  of  several  possible  types 
and  may  not  always  be  an  ideal  model,  it  is  still  relevant  to 
acquisition  inaccuracy  and  corruption  processes  that  may  be 
encountered  across  several  medical  imaging  modalities.  The 
use  of  an  intensity-based  similarity  metric  appears  to  give 
the  method  an  advantage  in  being  generally  insensitive  to 
reasonably  expected  amounts  of  image  noise. 

Figure  3  demonstrates  that  because  of  the  random 
nature  of  the  boundary  condition  errors,  the  magnitude  is 
itself  not  an  accurate  indicator  of  reconstruction  quality. 
This  necessitated  the  introduction  of  a  more  suitable 
parameter  that  accounts  for  the  net  effect  of  the  altered 
boundary  conditions  in  order  to  perform  fair  evaluations.  In 
essence,  randomizing  the  vectors  at  every  node  causes  the 
optimization  to  use  an  unpredictable  starting  pose  and 
increases  its  chance  of  converging  to  an  improper  minimum. 
Therefore,  the  ‘initial  alignment  error’  (AE)  is  defined  as  the 
relative  percent  change  between  the  objective  function 
evaluation  using  the  gold  standard  boundary  conditions  and 
those  of  the  test  case.  An  as  example,  it  could  be  assumed 
that  vectors  of  magnitude  0.5  would  be  a  much  more 
tolerable  error  than  2.0,  but  it  is  the  significantly  larger  AE 
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Figure  4.  Reconstruction  quality  vs.  percent  additive  image  noise. 
The  drop-off  after  10%  additive  noise  indicates  the  threshold  of 
tolerance  for  the  method. 


Figure  5.  Reconstruction  quality  vs.  percent  change  in  initial 
alignment  relative  to  gold  standard.  The  majority  of  errors  tested 
produced  satisfactory  reconstructions. 


of  the  former  that  actually  predicts  the  poor  outcome. 
However,  it  should  also  be  noted  (results  not  shown  here) 
that  even  if  the  same  set  of  error  vectors  are  scaled  over 
varying  magnitudes,  there  is  no  clear  trend  in  alignment 
error  or  reconstruction  quality.  This  appears  to  imply  that 
certain  boundary  nodes,  most  likely  those  in  the  direction  of 
largest  strain,  have  a  greater  effect  on  reconstruction  and 
merit  particular  care  in  selection.  Other  factors  influencing 
unfavorable  reconstructions  are  most  likely  nonlinear  effects 
not  predicted  by  the  current  model  as  well  as  an  inherent 
lack  of  discrimination  by  intensity-based  similarity  metrics 
in  analyzing  the  regularity  of  the  imposed  grid  pattern. 
Nevertheless,  for  the  error  magnitudes  tested  that  best 
approximate  realistic  inaccuracies  (i.e.  <0.5  units),  the 
alignment  errors  were  small  and  quality  of  the  end 
reconstruction  was  seen  to  be  quite  good.  This  qualitatively 
validates  the  current  method  of  determining  point 
correspondence  as  a  reasonable  procedure  with  an 
accommodating  margin  (factor  of  four)  in  light  of  typical 
user  interaction. 

5.  CONCLUSIONS 

In  this  work,  we  have  presented  a  method  for  recovering 
elasticity  parameters  from  image  data  of  thin  membrane 
structures  by  maximizing  the  image  similarity  between  two 
different  states  of  mechanical  loading  within  the  context  of 
an  inverse  problem.  The  biomechanical  model,  multi¬ 
resolution  optimization,  and  image  acquisition  are  each 
modular  components  of  this  elastographic  reconstruction 
framework,  making  it  extensible  to  added  sophistication. 
Tests  were  conducted  to  examine  the  tolerance  of  the 
method  to  degraded  or  improper  inputs.  The  results  indicate 
that  the  gold  standard  data  set  was  mostly  optimal  for 
obtaining  a  successful  reconstruction.  Widening  disparities 
in  either  image  data  or  boundary  condition  selection  from 
that  in  the  gold  standard  caused  observable  trends  of 
declining  reconstruction  quality.  Based  on  these 
observations,  it  appears  that  the  method  handles  most 


expected  variations  encountered  in  image  acquisition  as  well 
as  the  majority  of  typical  user  inaccuracies.  Because  there 
are  complicated  effects  associated  with  mapping  of  the 
Dirichlet  boundary  conditions  in  constraining  the 
displacement  solution  of  the  model,  further  study  on  inter¬ 
rater  variability  in  selection  as  well  as  comparisons  with 
more  automated  point  correspondence  methods  is  ongoing. 
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BACKGROUND 

RESULTS 

Changes  to  the  local  cytoarchitecture  induced  in  a  variety  of  pathologies  can  manifest 
as  alterations  in  tissue  elasticity  that  are  relevant  in  clinical  examination  and  evaluation. 
Many  elastography  methods  are  typically  dependent  on  the  specific  modality  around 
which  they  were  developed  (e.g.  magnetic  resonance  and  ultrasound  imaging).  We 
have  developed  ‘modality-independent  elastography’  (MIE)  as  a  reconstruction  method 
that  recovers  the  material  properties  of  soft  tissue  via  model-based  analysis  of  image 
data  acquired  at  different  states  of  mechanical  loading.  The  algorithm  utilizes  image 
similarity  in  the  performance  of  a  multi-resolution,  non-linear  optimization.  Previous 
work  with  a  phantom  membrane  used  for  simulated  dermoscopic  applications  prompted 
this  preliminary  investigation  of  the  relative  effects  of  additive  image  noise  and 
boundary  condition  determination  errors  on  the  performance  of  the  method.  The  results 
as  quantified  by  elasticity  contrast  and  localization  accuracy  indicate  that  the 
reconstruction  process  is  robust  in  the  presence  of  realistic  levels  of  image  corruption 
and  tolerates  the  majority  of  boundary  condition  mapping  errors. 


PURPOSE 


The  inputs  to  the  reconstruction  process  are  in  two  major  forms:  image  data  and 
boundary  condition  estimation.  Inadequate  fidelity  in  either  quantity  is  capable  of 
affecting  the  success  of  the  reconstruction  through  some  form  of  model-data  mismatch. 
We  proposed  to  test  the  sensitivity  of  the  algorithm  to  various  levels  of  an  applied  noise 
process  by  altering  either  the  intensity  distribution  of  the  target  image  or  the 
displacement  vectors  defining  the  Dirichlet  boundary  conditions. 


METHODS 


Figure  1 .  Flow  chart  of  MIE.  After  acquisition,  source  and  target  images  (A)  are  discretized  into  regions  and  zones,  respectively.  The  reconstruction  process 
involves  updating  elastic  modulus  values  (B,E)  to  drive  a  finite  element  model-based  image  deformation  (C)  until  the  best  match  is  found  (D). 
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Figure  2.  MIE  reconstruction  experiment. 

(Left  panel)  A  two-material  phantom  mimicking 
skin  was  constructed  as  a  thin  membrane 
measuring  15  cm  x  15  cm,  with  a  single  5  cm 
circular  stiff  inclusion  embedded  in  the  center.  The 
phantom  was  manufactured  with  like-colored 
polyurethanes  which  have  an  inclusion-to-bulk 
elasticity  contrast  of  approximately  5.7:1.  The 
membrane  was  stretched  in  a  uniaxial  fashion 
while  a  CCD  camera  mounted  above  acquired 
image  pairs  of  the  membrane  in  pre-  and  post- 
deformed  states  (960  x  1280  pixel  resolution,  8-bit 
grayscale). 

(Right  panel)  Top  row:  Source  and  target  images 
with  overiay  of  finite  element  mesh  boundaries 
(red)  that  demarcate  the  area  reconstructed. 
Below:  Reconstruction  progression  over  increasing 
number  of  regions  (N  =  16,  64,  256,  400)  to  refine 
the  spatial  distribution  of  elasticity  values.  This 
reconstruction  serves  as  the  gold  standard  for  the 
remainder  of  this  work. 


Figure  3.  MIE  image  noise  reconstruction  experiment. 

Gaussian  random  fields  of  variable  strength  with  respect  to  the 
variance  of  non-background  pixel  values  were  applied  in  an 
additive  fashion  to  the  target  image.  Shown  in  the  left  column 
from  top  to  bottom  are  the  original  target  and  then  with  10%, 
20%,  and  30%  noise.  In  the  right  column  are  the  corresponding 
elasticity  reconstructions  after  application  of  a  thresholding 
scheme  to  classify  bulk  (black)  and  inclusion  materials 
(white/gray).  The  known  segmentation  of  the  inclusion  was  used 
to  retrospectively  calculate  the  positive  predictive  value  of 
identifying  the  correct  material  type  within  the  proper  boundaries 
as  well  as  the  mean  elasticity  contrast  of  the  overall  distribution. 
For  this  work,  our  overall  evaluation  of  reconstruction  quality  is 
expressed  as  the  product  of  these  two  quantities.  The  effect  of 
additive  noise  is  to  decrease  reconstruction  quality  as  evidenced 
in  the  progressively  poorer  localization  of  the  inclusion. 


Figure  4.  MIE  boundary  condition 
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Randomized  vectors  of  a  particular  magnitude  were  applied  to  the  boundary  condition  specifications  of  the  same  finite  element  mesh  used  for  all  reconstructions.  This 
simulates  targeting  error  by  the  user  in  the  currently  semi-automated  method  of  point  correspondence  selection,  and  the  effect  is  illustrated  in  the  top  row:  from  left  to 
right,  the  gold  standard  boundary  and  then  with  mis-estimation  of  0.75,  1.0,  and  2.0  mesh  units  (equivalent  to  pixel  coordinates)  in  the  Dirichlet  conditions  (slightly 
exaggerated  scale  for  visual  effect).  The  corresponding  reconstructions  in  the  middle  and  bottom  rows  demonstrate  that  two  different  trials  using  the  same  magnitude  of 
randomized  vectors  can  effect  very  different  levels  of  reconstruction  quality. 


Figure  5.  Reconstruction  quality  vs.  image  noise. 

Three  trials  of  image  noise  were  performed  (shown  averaged  and  with 
standard  error  bars);  the  drop-off  in  reconstruction  quality  indicates 
the  presence  of  a  threshold  at  approximately  1 0%  additive  Gaussian 

Figure  6.  Reconstruction  quality  vs.  boundary  condition 

Two  trials  of  eight  levels  of  noise  ranging  from  0.1  to  2.0  mesh  units 
were  performed.  Each  reconstruction  was  treated  as  a  separate  data 
point  based  on  its  initial  alignment  error,  defined  here  as  the  relative 
change  between  the  objective  function  evaluation  using  the  gold 
standard  boundary  conditions  and  those  of  the  [randomized]  test 


Figures  3  and  5  show  that  the  achieved  reconstruction  becomes  more  inaccurate  with  increased  image  noise.  However,  the  ability  to 
identify  and  localize  the  stiff  inclusion  is  not  significantly  impaired  until  a  noise  field  of  greater  than  10%  is  applied.  The  threshold  was 
found  by  determining  which  level  of  noise  provided  the  best  minimum  sum  squared  error  fit  of  two  lines  to  the  distribution  of 
reconstruction  quality.  This  would  indicate  that  the  similarity  metric  and  objective  function  are  robust  to  intensity  deviations  of  about  6 
gray  levels  in  an  8-bit  image.  While  Gaussian  noise  may  not  always  be  an  ideal  model,  as  a  preliminary  point  of  investigation,  it  is  still 
relevant  to  acquisition  inaccuracy  and  corruption  processes  that  can  be  encountered  across  several  medical  imaging  modalities. 
Figure  4  demonstrates  that  the  magnitude  of  the  random  vectors  is  itself  not  an  accurate  indicator  of  reconstruction  quality  because 
the  multiple  degrees  of  freedom  afforded  by  the  boundary  nodes  cause  the  optimization  to  use  an  unpredictable  starting  pose, 
increasing  the  chances  of  converging  to  an  improper  local  minimum.  This  necessitated  the  introduction  of  the  initial  alignment  error 
(AE)  to  provide  a  consistent  means  of  comparison  between  trials  (Figure  6).  As  a  further  example,  it  could  be  assumed  that  vectors 
of  magnitude  0.5  would  be  a  much  more  tolerable  error  than  2.0,  but  it  is  the  significantly  larger  AE  of  the  former  that  actually 
predicts  the  poor  outcome.  It  should  also  be  noted  (results  not  shown  here)  that  even  if  the  same  set  of  error  vectors  are  scaled  over 
varying  magnitudes,  there  is  no  clear  trend  in  alignment  error  or  reconstruction  quality.  This  appears  to  imply  that  certain  boundary 
nodes,  most  likely  those  in  the  direction  of  largest  strain,  have  a  greater  effect  on  reconstruction  and  merit  particular  care  in 
selection.  Nevertheless,  for  the  error  magnitudes  that  approximate  inaccuracy  in  boundary  condition  demarcation  (i.e.  <0.5  units),  the 
quality  of  those  reconstructions  was  satisfactory.  This  qualitatively  validates  the  current  method  of  determining  point  correspondence 
as  a  reasonable  procedure  with  an  accommodating  margin  (factor  of  four)  in  light  of  typical  user  interaction.  Further  research  is 
ongoing  into  validation  and  control  of  boundary  conditions,  as  well  as  more  automated  methods  of  point  correspondence. 
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ABSTRACT 

A  significant  amount  of  breast  cancer  research  in  recent  years  has  been  devoted  to  novel  means  of  tumor  detection 
such  as  MR  contrast  enhancement,  electrical  impedance  tomography,  microwave  imaging,  and  elastography.  Many  of 
these  detection  methods  involve  deforming  the  breast.  Often,  these  deformed  images  need  to  be  correlated  to  anatomical 
images  of  the  breast  in  a  different  configuration.  In  the  case  of  our  elastography  framework,  a  series  of  comparisons 
between  the  pre-  and  post-deformed  images  needs  to  be  performed.  This  paper  presents  an  automatic  method  for 
determining  correspondence  between  images  of  a  pendant  breast  and  a  partially-constrained,  compressed  breast.  The 
algorithm  is  an  extension  to  the  symmetric  closest  point  approach  of  Papademetris  et  al.  However,  because  of  the  unique 
deformation  and  shape  change  of  a  partially-constrained,  compressed  breast,  the  algorithm  was  modified  through  the  use 
of  iterative  closest  point  (ICP)  registration  on  easily  identifiable  sections  of  the  breast  images  and  through  weighting  the 
symmetric  nearest  neighbor  correspondence.  The  algorithm  presented  in  this  paper  significantly  improves 
correspondence  determination  between  the  pre-  and  post-deformed  images  for  a  simulation  when  compared  to  the 
original  Papademetris  et  al.’s  symmetric  closest  point  criteria. 

Keywords:  deformation,  registration,  finite  element,  breast  cancer,  elastography,  soft-tissue  mechanics 

1.  INTRODUCTION 

In  a  2001  report  released  by  the  American  Cancer  Society,  a  woman  in  the  United  States  had  a  1  in  8  lifetime  risk  of 
developing  breast  cancer,  and  these  odds  increased  depending  on  factors  such  as  age,  family  history,  ethnicity,  and 
genetic  predisposition.  It  also  presented  statistics  illustrating  the  strong  connection  between  survival  rate  and  stage  of 
cancer  development  at  diagnosis  (localized  (96.4%),  regional  (77.7%),  and  distant  (21.1%))*[1].  Given  the  obvious 
benefits  of  extended  lifetime  with  early  diagnosis,  there  has  been  a  resultant  push  to  increase  the  sensitivity  and 
specificity  of  current  imaging  methods  and  modalities.  However,  in  addition  to  improving  upon  current  tumor-detection 
methodologies,  much  recent  research  has  been  devoted  to  novel  imaging  techniques  devised  to  supplement  diagnosis. 

Many  of  these  new  techniques  are  based  upon  the  exploitation  of  constitutive  properties  of  tissue,  e.g.  electrical, 
optical,  or  stiffness  properties.  The  differences  between  healthy  and  cancerous  tissue  with  regards  to  these  properties  has 
been  investigated  to  varying  degrees  for  a  number  of  years  and  suggested  results  have  been  forthcoming  [2-4]. 
Currently,  the  use  of  the  differences  in  these  properties  as  a  tumor  detection  method  can  be  seen  through  microwave 
imaging  [5],  electrical  impedance  tomography  (EIT)  [6],  near- infrared  tomography  [3],  ultrasound  elastography  (USE), 
magnetic  resonance  elastography  (MRE),  and,  most  recently,  in  modality  independent  elastography  (MIE)  [7].  Many  of 
these  detection  methods,  particularly  elastography,  involve  deforming  the  breast  and  comparing  the  pre-  and  post- 
deformed  images  of  the  breast.  In  addition,  there  are  other  methods  whereby  deformations  occur  due  to  repositioning  of 
the  patient  as  in  MR  contrast  enhancement  [8]. 

The  purpose  of  this  paper  is  to  present  an  algorithm  that  automatically  determines  correspondence  between  images 
of  the  pendant  and  deformed  breast  geometry.  This  algorithm  is  based  on  the  symmetric  nearest  neighbor  algorithm 
presented  by  Papademetris  et  al.  [9];  however,  the  algorithm  presented  in  this  paper  uses  certain  geometric  characteristics 
of  the  breast  to  improve  upon  the  accuracy  of  Papademetris  et  al.’s  method.  Due  to  the  unique  deformation  and  shape 
These  statistics  represent  the  5-year  survival  rate. 
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change  of  a  partially-constrained,  compressed  breast,  the  algorithm  was  modified  through  the  use  of  iterative  closest 
point  (ICP)  registration  on  easily  identifiable  sections  of  the  breast  images,  and  it  included  weighting  the  symmetric 
nearest  neighbor  correspondence.  The  work  presented  here  has  implications  within  the  modality  independent 
elastography  method  presented  by  Miga  et  al.  [7,  10]  as  well  as  for  the  multi-modal  registration  among  alternative 
imaging  techniques  for  breast  cancer. 

2.  MATERIALS  AND  METHODS 

2.1.  Breast  surface  correspondence  algorithm 

A  model  of  a  human  breast  is  generated  in  the  pendant  position  based  on  breast  CT  data  (courtesy  of  John  M. 
Boone,  Ph.D.,  Professor  of  Biomedical  Engineering,  University  of  California  at  Davis  [11-12]).  Using  a  finite  element 
model,  a  deformation  simulating  the  inflation  of  a  membrane  pressing  against  the  breast  is  performed,  creating  a 
deformed  breast  geometry.  Correspondence  between  point  clouds  of  the  pendant  breast  and  the  compressed  breast  is 
performed  using  a  multi-step  procedure  based  on  the  symmetric  closest  point  algorithm  of  Papademetris  et  al.  The 
algorithm  was  modified  such  that  easily  identifiable  regions  are  constrained  using  an  iterative  closest  point  (ICP) 
registration  to  initially  establish  correspondence.  End  correspondence  of  a  point  is  still  partially  determined  by  its 
neighbors’  correspondence.  Thus,  the  algorithm  consists  of  four  main  steps:  ICP  registration  of  the  base  and  nipple 
regions  of  the  breast,  the  initial  correspondence  determination  by  symmetric  nearest  neighbors,  the  “propagating  front” 
of  the  correspondence  through  the  remaining  points  in  the  point  cloud,  and  the  final  smoothing  and  mapping  of  the 
correspondence. 

For  purposes  of  simplicity,  the  following  nomenclature  will  be  used  in  describing  the  algorithm:  the  point  cloud  that 
defines  the  pendant  breast  is  the  source  point  cloud,  s2,  and  the  point  cloud  that  defines  the  partially-constrained, 
compressed  breast  is  the  target  point  cloud,  s2.  Also,  a  “faux  registration”  of  sj  to  s2  means  that  the  rotation  and 
translation  matrices  necessary  to  ICP  register  sj  to  s2  are  known,  but  the  actual  registration  of  sj  to  s2  does  not  take  place. 
Furthermore,  the  algorithm  is  constructed  for  a  super-sampled  target  point  cloud. 

2.1.1.  Initial  correspondence  via  ICP  partial  volume  registration 

The  purpose  of  this  step  is  to  use  the  more  prominent  geometric  characteristics  of  the  breast  and  the  fixation  of  the 
breast  at  the  chest  wall,  namely  the  regions  around  the  nipple  and  the  base  of  the  breast,  as  a  location  for  accurately 
determining  initial  correspondence.  However,  due  to  the  compression,  the  nipple  of  the  deformed  breast  (the  target  point 
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Figure  1 :  Point  clouds  of  the  deformed  and  pendant  breast,  with  the  top  20%  of  the  pendant  breast  registered  to  the  deformed  breast 
according  to  the  rotation  and  translation  matrices  provided  from  a  faux  registration  of  the  top  46.25%  of  the  breasts. 
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cloud)  is  translated  and  lies  at  a  different  angle  than  the  nipple  of  the  pendant  breast  (the  source  point  cloud).  By 
minimizing  the  average  distance  from  the  points  in  the  nipple  region  on  the  registered  source  point  cloud  to  the  points  in 
the  nipple  region  on  the  target  point  cloud,  the  nipple  regions  of  the  pendant  and  deformed  breast  are  most  accurately 
aligned.  (The  average  distance  means  the  average  Euclidean  distance  between  a  point  on  the  source  point  cloud  and  its 
closest  point  on  the  target  point  cloud.)  For  the  nipple  region  registration  (which  the  nipple  region  is  assumed  to  reside 
at  the  upper  20%  of  the  breast’s  height),  there  exists  a  certain  height  percentage,  y/,  such  that  if  all  points  of  the  pendant 
breast  at  or  above  yj  were  faux  registered  to  all  of  the  points  of  the  deformed  breast  at  or  above  y/,  then  the  resultant 
rotation  and  translation  matrices  would  minimize  the  distance  between  nipple  regions  of  the  pendant  and  deformed 
breast.  This  height  percentage,  iff,  is  determined  as  follows:  a  binary  search  algorithm  faux  registers  a  range  of  height 
percentages,  from  the  upper  50%  to  the  upper  20%  of  the  breast  height.  The  algorithm  then  compares  the  average 
distance  between  the  nipple  regions  for  the  rotation  and  translation  matrices  produced  by  the  faux  registration.  The 
height  percentage  that  produces  the  minimal  distance  between  nipple  regions  is  designated  y/.  Finally,  by  using  the 
rotation  and  translation  matrices  defined  by  the  registration  of  the  y/  percentage,  the  nipple  region  of  the  source  point 
cloud  is  registered  to  the  nipple  region  of  the  target  point  cloud  (see  Figure  1),  and  initial  correspondence  between  the 
point  clouds  can  now  take  place.  This  process  is  repeated  for  the  base  regions  (i.e.  the  lower  5%  of  the  breast  height)  of 
the  source  and  target  point  clouds. 

2.1.2.  Symmetric  nearest  neighbor  correspondence 

Initial  correspondence  is  restricted  to  the  registered  regions  of  the  point  clouds  (the  nipple  and  base  regions). 
Correspondence  is  achieved  through  a  method  Papademetris  et  al.  coined  “symmetric  nearest  neighbor.”  For  a  point  p3 
on  the  source  point  cloud  sj,  the  closest  point  to  pj  on  the  target  point  cloud  s2  is  p2.  In  other  words,  p2  is  chosen  such 
that  the  distance  from  p3  to  p2  is  the  shortest  Euclidean  distance  possible  when  compared  to  the  Euclidean  distance 
between  p3  and  any  other  point  on  s2.  p2  then  determines  its  closest  point  on  s3  to  be  p.  If  p2=  p,  then  the  points  p3  and 
p2  are  corresponding  and  called  symmetric  nearest  neighbors.  The  vector  from p3  to  p2  is  called  a  corresponding  vector. 


Figure  2:  Simple  diagram  of  a  “propagating  front”  of  correspondence.  The  weighted  average  of  the  symmetric  closest  point  vectors 
from  C],  c2,  c3  to  bh  b2 ,  b3  create  vAVE.  p3  is  a  point  that  does  not  have  correspondence,  p  is  created  by  projecting  along  vave  from ph 
and  the  closest  point  p2  to  p  becomes  p{  s  corresponding  point  (with  corresponding  vector  vCOR1 ?).  Correspondence  then  bleeds  to  T, 
which  now  uses  vCORR  in  the  creation  of  its  weighted  average  vector. 


2.1.3.  “Propagating  front”  of  correspondence 

Correspondence  then  continues  to  be  established  using  a  “propagating  front”  methodology  in  which  points  on  the 
source  point  cloud  that  have  neighbors  with  a  high  degree  of  correspondence  are  the  first  candidates  for  establishing 
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correspondence,  as  shown  in  Figure  2.  For  point  /?7  on  sj,  the  first-  and  second-order  neighboring  points  of  pj  in  Sj  that 
already  have  correspondence  are  c7,  c2,...cn  (points  bh  b2,...bn  are  the  corresponding  points  of  c7,  c2,...cn  on  the  target 
point  cloud  s2,.and  their  corresponding  vectors  are  v7,  v2,...vn).  An  average  of  the  corresponding  vectors  v7,  v2,...vm 
weighted  by  the  distance  of  c7,  c2,...cn  to  /?7,  is  then  computed  and  is  denoted  vAVE.  A  temporary  point,  p,  is  then 
projected  by  vAVE  from  /?7,  and  the  closest  point  (as  determined  by  Euclidean  distance)  on  s2  to  p  is  p2.  p2  is  the 
corresponding  point  for  /?7.  The  vector  from p!top2  is  the  corresponding  vector,  vCORR.  p2  can  be  any  point  in  s2  unless  it 
is  bh  b2,...bn.  In  this  case,  the  next  closest  point  to  P  on  ^  that  is  not  b^  b2,...bn  is  the  corresponding  point  to />7.  This 
step  will  provide  correspondence  for  every  point  in  the  source  point  cloud  as  long  as  there  is  at  least  one  initial  point 
correspondence  between  the  source  and  target  point  clouds.  The  algorithm  continues  in  an  iterative  fashion  such  that 
correspondence  will  be  determined  only  for  those  points  in  sj  with  a  maximum  number  of  neighbors  with 
correspondence. 

2.1.4.  Smoothing  and  mapping 

After  every  point  in  the  source  point  cloud  has  correspondence,  the  entire  point  cloud  undergoes  a  smoothing  and 
mapping  to  reach  a  final  correspondence.  Point  /?7  of  source  point  cloud  sj  has  neighboring  points  c7,  c2 ,. .  ,c„  .  Points  c7, 
c2,...cn  have  corresponding  points  bh  b2,...bn  in  the  target  point  cloud  s2 ,  and  their  corresponding  vectors  are  v7,  v2,...vn. 
An  average  of  the  corresponding  vectors  v7,  v2,...vn,  weighted  by  the  distance  of  c7,  c2,...cn  to  ph  is  computed  and 
denoted  vAVE.  A  temporary  point,  P,  is  then  projected  by  vAVE  from  p 7,  and  the  closest  point  (as  determined  by  Euclidean 
distance)  from  s2  to  P  is  p2  ,  which  is  the  final  corresponding  point  of  pj. 

2.2.  Breast  deformation  simulation 

In  order  to  test  the  algorithm  described  above,  a  finite  element  (FE)  model  of  a  breast  was  generated  based  on  a 
subject’s  CT  image  volume.  The  biomechanics  of  breast  deformation  was  described  using  a  three-dimensional  linear 
elastic  Hookean  solid.  Once  the  geometry  was  generated,  boundary  conditions  were  implemented  that  simulated  the 
compression  of  the  breast  by  an  adjacent  balloon  inflation  cuff.  Figure  3  shows  the  geometry  and  the  boundary 
conditions  applied.  Results  from  varying  the  strength  of  the  stress  field  can  be  seen  in  Figure  4.  Figure  4a  shows  the 
pre-post  compression  of  the  breast  subject  to  the  boundary  conditions  in  Figure  3.  Figure  4b  and  4c  represent  different 
inflation  states,  75%  and  50%  respectively.  By  using  a  finite  element  model,  exact  correspondence  between  source  and 
target  can  be  established  and  used  to  rate  the  algorithm  performance.  It  should  be  noted  though  that  the  target  cloud,  i.e. 
deformed  finite  element  breast  mesh,  was  super-sampled  to  allow  for  more  variability  when  executing  the 
correspondence  algorithm.  In  the  clinical  context,  the  point-cloud  representations  would  represent  two  unique  image 
acquisitions  so  a  corresponding  structured  grid  (as  provide  by  the  FE  mesh)  would  not  exist. 
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Figure  3.  Finite  element  breast  model  with  boundary  conditions  shown.  Grayscale  represents  the  application  of  a  normal  stress 
in  kPa.  The  base  of  the  volume  was  fixed. 
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Figure  4:  Pendant  and  deformed  breast  point  clouds  for  a.  100%  simulation  b.  70%  simulation  and  c.  50%  simulation. 


3.  RESULTS 
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Figure  5.  (a)  Distribution  of  correspondence  error  using  the  proposed  algorithm  with  the  100%  deformation  FE  simulation, 
(b)  Distribution  of  correspondence  error  with  the  70%  simulation,  (c)  Distribution  of  correspondence  error  with  the  50% 
simulation,  (d)  Distribution  of  correspondence  error  for  100%  deformation  FE  simulation  using  Papademetris  et  al.’s  algorithm. 
Grayscale  shown  is  common  to  all  figures. 


29 


Figure  5a-c  shows  the  result  of  running  the  various  simulations  through  our  algorithm  while  Figure  5d  shows  the 
result  of  running  the  100%  simulation  with  using  solely  the  symmetric  closest  point  criteria  of  Papademetris  et  al.’s 
algorithm.  The  shading  of  each  deformed  breast  is  root  mean  square  (RMS)  error  of  the  distance  from  the  predicted 
target  corresponding  point  and  the  actual  target  corresponding  point.  Figure  6  compares  the  RMS  error  and  maximum 
error  (again,  the  distance  from  the  predicted  target  corresponding  point  to  the  actual  target  corresponding  point)  of  our 
algorithm  versus  the  maximum  displacement  compression  exerted  on  the  breast  for  each  simulation.  Additionally,  the 
maximum  error  for  Papademetris  et  al.’s  algorithm  and  the  100%  simulation  is  2.78  cm,  and  the  RMS  error  is  1.02  cm. 


Maximum  compression  of  breast  (meters) 


Figure  6:  Comparison  of  the  maximum  displacement  in  the  depression  of  the  breast  to  the  RMS  error  (•)  and  the  maximum  error 
(■)  for  the  100%  simulation  (right  points),  70%  simulation  (middle  points)  and  50%  simulation  (left  points). 


4.  DISCUSSION 

Figure  5a-c  demonstrates  that  our  algorithm,  overall,  assigns  correspondence  from  the  pendant  to  the  deformed 
breast  reasonably  well,  especially  at  less  severe  deformations.  Figure  6  provides  further  support,  demonstrating  that  the 
relatively  low  RMS  error  ranges  from  1.2  mm  for  the  50%  simulation  to  6.0  mm  for  the  100%  simulation.  Furthermore, 
there  is  a  linear  relationship  between  degree  of  compression  and  the  maximum  and  RMS  error  for  our  algorithm,  with 
both  errors  decreasing  as  the  degree  of  compression  decreases.  However,  there  is  still  error  in  the  correspondence 
assignments,  primarily  at  the  site  of  compression.  This  error  is  due  to  the  algorithm’s  inability  to  cope  with  the  steep 
point  translations  associated  with  the  depression;  the  algorithm  continues  to  assign  correspondence  to  the  edge  of  the 
depression  rather  than  penetrating  into  the  depression.  This  suggests  that  further  research  into  the  incorporation  of  some 
other  type  of  information  into  our  algorithm  is  needed  for  greater  accuracy.  But,  as  evidenced  by  Figure  5d,  our 
algorithm  is  significantly  more  accurate  than  the  original  Papademetris  et  al.  algorithm.  While  the  Papademetris  et  al. 
algorithm  handles  the  depression  approximately  as  well  as  our  algorithm  with  a  maximum  error  of  2.78  cm  (compared  to 
our  algorithm’s  maximum  error  of  2.68  cm  for  the  same  simulation),  the  accuracy  of  the  correspondence  assignments 
continues  to  decrease  up  into  the  nipple  region.  This  results  in  a  RMS  error  of  1.02  cm,  which  is  significantly  higher 
than  our  RMS  error  of  6.0  mm.  It  should  be  noted  though  that  only  the  symmetric  closest  point  portion  of  the 
Papademtris  et  al.  algorithm  has  been  implemented.  The  algorithm  in  its  entirety  also  utilizes  a  curvature  mapping 
aspect  for  their  particular  application  in  cardiac  deformation  mapping.  Usage  of  curvature  information  would  not  be 
useful  for  breast  applications  due  to  its  geometrically  homogeneous  nature. 
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The  main  reason  our  algorithm  is  more  accurate  than  just  utilizing  the  symmetric  closest  point  algorithm  is  due  to 
the  first  step,  the  initial  registration  of  the  nipple  and  base  regions,  which  constrains  the  search  greatly.  By  rigidly 
registering  the  nipple  and  base  regions,  the  initial  pass  of  establishing  correspondence  is  highly  accurate.  In  the  original 
algorithm,  the  initial  pass  of  determining  correspondence  frequently  was  inaccurate  (the  entire  point  cloud  was  permitted 
to  initially  establish  correspondence),  particularly  at  the  nipple  region.  Because  the  deformed  breast’s  nipple  was 
translated  away  from  the  pendant  breast’s  nipple,  correspondence  was  established  between  the  far  side  of  the  pendant 
breast’s  nipple  and  the  near  side  of  the  deformed  breast’s  nipple  as  opposed  to  the  far  side  of  the  deformed  breast’s 
nipple. 

At  this  time,  it  needs  to  be  noted  that  this  paper’s  algorithm  needs  further  study  with  respect  to  the  number  and  type 
of  test  cases  (this  algorithm  assumed  the  nipple  resided  in  the  upper  20%  of  the  breast’s  height,  which  may  not  always  be 
the  case).  The  ICP  registration  of  the  nipple  region  may  introduce  an  error  when  our  algorithm  is  used  with  other 
simulations.  ICP  registration  may  produce  a  rotation  causing  the  source  point  cloud  to  spin  around  the  longitudinal  axis 
of  the  breast  (a  full  180°  in  the  most  extreme  case)  resulting  in  inaccurate  correspondence.  Given  the  intention  of  the 
step  is  to  provide  accurate  correspondence  of  the  nipple  (an  easily  identifiable  structure),  a  more  robust  algorithm  that 
prohibits  significant  rotation  may  be  needed.  In  addition,  we  should  also  note  that  without  careful  validation  of  the  FE 
model  for  deformation,  it  may  be  that  the  boundary  conditions  represented  by  Figure  3  do  not  accurately  describe  the 
displacements  associated  with  an  inflating  adjacent  balloon. 

5.  CONCLUSIONS 

The  addition  of  rigid  registration  initialization  of  the  chest-wall  and  nipple  region  of  the  breast  to  the  Papademetris 
et  al.  symmetric  closest  point  criteria  significantly  improves  the  accuracy  of  the  approach  for  the  breast’s  morphology. 
While  the  maximum  error  located  in  the  center  of  the  depression  remains  high  for  our  algorithm  at  2.68  cm  (similar  to 
Papademetris  et  al.’s  maximum  error  of  2.78  cm),  the  RMS  error  is  significantly  improved  to  6.0  mm  (as  compared  to 
Papademetris  et  al.’s  RMS  error  of  1.02  cm).  Furthermore,  it  is  clear  that  the  algorithm  demonstrates  an  increase  in 
accuracy  as  the  degree  of  compression  is  decreased.  Finally,  the  use  of  other  geometric  information  and  additional 
experimentation  with  test  cases  may  enhance  this  algorithm  to  provide  more  robust  and  accurate  determination  of 
correspondence. 
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ABSTRACT 

Multiple  skin  conditions  exist  which  involve  clinically  significant  changes  in  elastic  properties. 
Early  detection  of  such  changes  may  prove  critical  in  formulating  a  proper  treatment  plan.  However, 
most  diagnoses  still  rely  primarily  on  visual  inspection  followed  by  biopsy  for  histological  analysis.  As  a 
result,  there  would  be  considerable  clinical  benefit  if  a  noninvasive  technology  to  study  the  skin  were 
available.  The  primary  hypothesis  of  this  work  is  that  skin  elasticity  may  serve  as  an  important  method 
for  assisting  diagnosis  and  treatment.  Perhaps  the  most  apparent  application  would  be  for  the 
differentiation  of  skin  cancers,  which  are  a  growing  health  concern  in  the  United  States  as  total  annual 
cases  are  now  being  reported  in  the  millions  by  the  American  Cancer  Society.  In  this  paper,  we  use  our 
novel  modality  independent  elastography  (MIE)  method  to  perform  dermoscopic  skin  elasticity 
evaluation.  The  framework  involves  applying  a  lateral  stretching  to  the  skin  in  which  dermoscopic 
images  are  acquired  before  and  after  mechanical  excitation.  Once  collected,  an  iterative  elastographic 
reconstmction  method  is  used  to  generate  images  of  tissue  elastic  properties  and  is  based  on  a  two- 
dimensional  (2-D)  membrane  model  framework.  Simulation  studies  are  performed  that  show  the  effects 
of  three-dimensional  data,  varying  subdermal  tissue  thickness,  and  nonlinear  large  deformations  on  the 
framework.  In  addition,  a  preliminary  in  vivo  reconstruction  is  demonstrated.  The  results  are 
encouraging  and  indicate  good  localization  with  satisfactory  degrees  of  elastic  contrast  resolution. 

KEY  WORDS:  Elastography,  Dermoscopy,  Mechanical  Properties,  Finite  Element,  Image  Similarity, 
Elasticity  Imaging 


1.  INTRODUCTION 

Because  multiple  conditions  exist  which  involve  clinically  significant  changes  in  skin  elasticity, 
early  detection  of  such  changes  could  prove  critical  in  formulating  a  proper  treatment  plan.  However, 
most  diagnoses  still  rely  primarily  on  visual  inspection  followed  by  biopsy  of  suspect  areas  for 
histological  analysis.  Perhaps  the  most  apparent  application  would  be  for  the  differentiation  of  skin 
cancers,  which  are  a  growing  health  concern  in  the  United  States  as  total  annual  cases  are  now  being 
reported  in  the  millions  by  the  American  Cancer  Society1.  Of  the  three  major  types  of  skin  cancer 
reported  annually,  basal  cell  carcinoma  (BCC)  makes  up  approximately  800,000+,  squamous  cell 
carcinoma  (SCC)  cases  number  approximately  200,000+,  and  a  remaining  60,000+  melanoma  cases.  With 
respect  to  BCC,  approximately  5-10%  of  these  can  be  aggressive  and  infiltrate  the  surrounding  tissue  and 
sometimes  into  bone  and  cartilage.  It  rarely  metastasizes  but  can  cause  scars  and  disfigurement.  With 
respect  to  SCC,  early  detection  is  the  key  to  successful  treatment.  If  left  unchecked,  SCC  can  also  cause 
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disfigurement  and  typically  approximately  3-4%  of  cases  results  in  metastasis  which  is  usually  fatal. 
Melanoma  is  the  most  fatal.  Melanoma  is  malignant  and  if  left  unchecked,  it  will  spread  to  other  parts  of 
the  body,  becomes  difficult  to  treat,  and  can  be  fatal.  If  one  were  to  include  among  the  common  cancer 
statistics,  aggressive  BCC’s  and  metastatic  SCC’s,  skin  cancer  would  likely  be  the  second  most  prevalent 
among  newly  diagnosed  cancers. 

While  the  above  statistics  are  compelling,  when  speaking  to  lethality,  skin  cancer  is  less 
significant  than  its  more  fatal  counterparts  in  breast,  lung,  colon/rectal,  and  prostate.  However  when 
considering  the  economic  cost  that  skin  cancer  incurs  on  US  healthcare,  the  pursuit  of  skin  cancer 
characterization  has  considerable  merit.  The  ability  to  differentiate  benign  from  malignant,  and 
aggressive  from  non-aggresssive  skin  lesions  would  provide  considerable  savings  to  health  care  costs 
within  the  United  States.  Each  year  it  is  estimated  that  approximately  5-7  million  patients  undergo 
biopsies  for  skin  cancer  with  only  a  fraction  resulting  in  malignancies.  While  complete  multi-center 
biopsy  studies  have  not  been  performed,  one  study  that  took  place  documented  the  percentage  of  skin 
biopsy  specimens  that  were  malignant,  i.e.  they  termed  this  the  “malignancy  ratio”  2.  In  this  study,  the 
malignancy  was  approximately  40%  with  a  very  wide  variability  among  the  22  dermatologists  (13.4- 
86.6%)  that  participated  in  the  study.  If  one  considers  5  million  biopsies,  this  would  translate  to 
approximately  3  million  unnecessary  biopsies  per  year,  and  approximately  $300-600  million  in  saved 
expenditures  per  year  3’ 4.  Furthermore,  this  does  not  even  include  cosmetic  surgery  costs  in  the  case  of 
scarring  or  complications  associated  with  bleeding  and  infections.  If  an  inexpensive  imaging  device 
could  differentiate  lesions  with  reasonably  high  specificity  and  sensitivity,  it  would  have  considerable 
significance. 

It  should  also  be  noted  that  measuring  cutaneous  elasticity  is  also  potentially  important  to  medical 
areas  outside  of  clinical  dermatology.  For  instance,  in  a  recent  study  of  100  women  receiving  hormone 
replacement  therapy,  Pierard  et  al  demonstrated  a  positive  correlation  between  bone  mass  density  and 
skin  elasticity  5’  6.  Another  study  conducted  by  Yoon  et  al  showed  a  similar  relevance  for  patients 
afflicted  with  diabetes  mellitus  7.  Further  uses  for  evaluating  skin  elasticity  range  from  surgery 
(minimization  of  scarring,  chronic  graft  versus  host  disease)  to  rheumatology  (scleroderma,  lupus)  to 
obstetrics  (striae  development  in  pregnancy) 8'12. 

With  respect  to  diagnostic  technological  advances,  developments  have  been  concerned  with 
obtaining  a  better  view  of  the  skin,  either  via  improved  optics  (i.e.  dermoscopy)  or  by  more  advanced  and 
novel  imaging  systems  ranging  from  high-frequency  ultrasound  to  confocal  laser  microscopy  13, 14.  Other 
strategies  involving  electrical  impedance  mismatch  15,  Raman  spectroscopy  16,  and  cytological  smears  14 
have  also  been  forthcoming.  In  lieu  of  these  methods  that  capitalize  on  electrical,  optical,  and 
biochemical  phenomena,  we  have  chosen  to  pursue  an  approach  which  is  based  on  analyzing  mechanical 
behavior  of  the  skin.  Detecting  changes  in  tissue  by  palpation  and  then  associating  them  with  a  disease 
state  has  a  long-standing  history  in  clinical  medicine,  and  utilizing  changes  in  the  mechanical  properties 
to  specifically  characterize  the  skin  does  have  precedent.  A  thoughtful  review  by  Edwards  and  Marks 
discusses  the  complex  mechanical  behavior  of  skin  when  subjected  to  in  vitro  and  in  vivo  testing  17.  Their 
review  highlights  extensive  methodologies  being  used  to  quantify  skin  properties  (e.g.  uni-axial  and  bi¬ 
axial  extensometry,  torsion  stimulators,  indentometery,  ballistometric  tests,  shear  wave  application 
devices,  dynamic  suction  methods,  ultrasonics,  and  electrodynamometry)  and  also  emphasizes  the 
difficulties  in  comparing  across  these  methods.  One  of  the  authors'  conclusions  is  that  the  need  for 
quantitative  reproducible  methods  to  assess  skin  health  is  necessary  given  the  considerable  subjectivity  in 
clinical  analysis. 

Following  previous  work  in  18,  we  are  developing  a  new  method  we  have  termed  ‘  modality  - 
independent  elastography’  (MIE)  that  combines  image  processing  with  an  inverse  problem.  More 
specifically,  image  similarity  metrics  routinely  used  with  image  registration  methods  are  recast  to  satisfy 
an  inverse  elasticity  problem  framework  whereby  mechanical  properties  within  a  biomechanical  model  of 
deforming  tissue  become  the  driving  parameters  for  improved  image  similarity.  In  this  way,  MIE 
circumvents  two  potential  limitations  of  current  elastographic  techniques.  First,  it  is  not  inherently 
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dependent  on  pre-processing  steps  such  as  homologous  feature  selection  and  tracking  which  drive  active 
contour  models  '  and  other  traditional  displacement-based  iterative  methods  '  ,  although  such 
techniques  can  aid  in  the  determination  of  boundary  conditions.  Secondly,  because  it  is  an  image 
processing  methodology,  MIE  is  not  reliant  on  a  particular  imaging  modality  (such  as  in  ultrasound  and 
magnetic  resonance  elastography)  as  long  as  the  acquired  images  provide  sufficient  pattern  to  allow  for 
comparison.  Building  upon  recently  with  a  multi-resolution  implementation25,  this  paper  presents 
analysis  using  a  tissue  model  that  incorporates  geometric  nonlinearities,  the  effects  of  the  three- 
dimensional  nature  of  the  problem  which  include  varying  subcutaneous  layer  thicknesses,  and  varying 
lesion  depth.  In  addition,  a  relatively  crude  preliminary  in  vivo  result  is  also  demonstrated. 

2.  MATERIALS  AND  METHODS 
2.1.  Modality  Independent  Elastography 


Fig  1.  MIE  framework. 

MIE  begins  with  the  acquisition  of  an  image  series  consisting  of  an  image  acquired  prior  to  and 
after  an  applied  deformation.  The  method  is  “independent”  because  it  does  not  require  any  specific 
imaging  modality  but  rather  that  sufficient  pattern  is  present  within  the  acquired  images  such  that  material 
properties  can  be  assessed  from  pattern  changes  within  the  acquisition-pair.  This  is  a  similar  requirement 
for  many  methods  of  non-rigid  image  registration.  At  its  core  MIE  is  an  image  analysis  method  whereby 
a  model-generated  deformation  field  is  applied  to  the  pre-deformed  image  series  {source)  and  compared 
to  its  acquired  deformed  counterpart  {target).  This  comparison  is  nested  within  a  Levenberg-Marquardt 
optimization  routine  such  that  the  material  properties  become  the  parameters  of  interest  in  matching  the 
model-deformed  source  image  to  the  acquired  target  image.  The  methods  and  development  of  this 
technique  have  been  reported  in  detail  elsewhere18,25'29. 

With  respect  to  the  optimization  framework  for  MIE,  it  can  be  represented  as  a  least  squared  error 
objective  function: 

4>(e)=  min{|s(E,)- s(ee|2}  (1) 

where  s(Et )  is  the  similarity  value  achieved  when  comparing  the  target  image  to  itself  (i.e.  the  maximum 

value  for  the  similarity  measure,  so  for  the  correlation  coefficient  S^Et)=  1)  and  s(ee)  is  the  similarity 
between  the  model-deformed  source  image  and  the  acquired  target  image  using  the  current  estimate  of  the 
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elastic  modulus.  Eq.  (1)  is  optimized  by  setting  the  partial  derivative  equal  to  zero  and  solved  using  a 
Levenberg-Marquardt  approach: 


[[JT][J]  +  a[l]]^E}=  [jt^(§,)-S(ee)}  (2) 

dSiE 

is  the  M  x  N  Jacobian  matrix  of  the  form  J  = - W 

8E 

measurements  (i.e.  zones)  being  made  and  N  is  the  number  of  material  property  regions.  Because 
NT][J]  (an  approximation  to  the  Hessian  matrix)  tends  to  be  ill-conditioned,  it  is  regularized  with  an 
empirically  determined  a  parameter  found  in  the  standard  Levenberg-Marquardt  approach30.  The 
determination  of  this  regularization  parameter  is  described  in  31.  The  multi-resolution  framework  to  MIE 
is  shown  in  Fig.  1.  The  methodology  used  involves  a  hierarchical  material  parameter  resolution  series. 
This  has  been  reported  elsewhere  25, 27  and  has  been  shown  to  assist  in  avoiding  local  minima  that  are 
associated  with  the  decorrelation  of  patterned  data. 

2.2.  Model  for  Skin 


where  [J] 


^  and  M  is  the  number  of  similarity 


One  critical  component  within  all  model-based  inverse  problem  frameworks  is  the  selection  of 
the  computer  model  to  represent  the  continuum  of  interest.  In  previous  phantom,  simulation,  and  in  vivo 
studies,  we  have  elected  to  employ  a  plane  stress  linear  elastic  model  to  simulate  the  skin.  This  model  is  a 
two  dimensional  approximation  of  the  three  dimensional  system  which  assumes  a  symmetric,  isotropic, 
thin  specimen  in  equilibrium,  and  stresses  that  are  constrained  to  lie  within  the  plane.  These  assumptions 
simplify  Cauchy’s  law  from  36  stiffness  constants  to  2  and  use  the  equation: 

V  •  a  =  0  (3) 

where  a  is  the  2D  Cartesian  stress  tensor  and  is  defined  below  followed  by  the  plane  stress  constitutive 
equations: 


u 

i— 

1  V  0 

au/dx 

b 

v  1  0 

dv/3 y 

ay 

"(1-v)2 

T 

o  o  1“V 

2 

dv/dx  +  du/dy 

where  E  is  Young’s  modulus,  and  v  is  Poisson’s  ratio.  For  this  work,  we  have  determined  that  a  Poisson’s 
ratio  of  0.485  for  our  skin  phantoms  and  tissue  work  has  performed  reasonably  well.  This  value  would 
correlate  with  a  -32:1  ratio  of  the  Lame'  constants  (k  :p,  with  pthe  shear  modulus,  and  X  being  the 
second  Lame'  constant)  which  is  reasonably  below  the  convention  for  Poisson  locking  (sometimes  called 
mesh  locking  and  typically  has  k  »  \i )  although  one  could  argue  that  hyperelastic  models  may  be  the 
more  appropriate  model  and  will  need  to  be  explored  in  the  future.  Using  the  Galerkin  weighted  residual 
method  to  integrate  and  solve  this  set  of  partial  differential  equations,  a  finite  element  framework  is 
generated  and  can  be  solved  to  represent  a  displacement  field  for  a  given  distribution  of  Young’s 
modulus. 

In  this  paper,  the  plane  stress  model  has  been  enhanced  to  incorporate  the  full  Green-Lagrange 
strain  tensor  as  defined  by: 


s= :  = 


Su,  dUj  Suk  3uk 


■  +  ■ 

dx,  dx ,  dx,  dx. 


(6) 


J  J 


Use  of  this  tensor  description  abandons  the  traditional  small-strain  approximation  in  favor  of  one 
compatible  with  large  deformations.  The  difference  in  solutions  between  small  and  large  deformation 
theory  can  be  seen  in  the  2D  simulations  shown  in  Fig.  2.  Fig.  2  compares  the  boundary  shape  of  a 
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compressive  strain. 

Hookean  linear  elastic  membrane  to  a  2D  Hookean  geometrically  nonlinear  elastic  model.  For  these 
simulations,  a  plane  stress  approximation  was  performed  for  comparison  only  of  the  solution  due  to 
nonlinear  terms  (i.e.  a  -30-40%  compressive  stress  could  never  be  applied  to  a  thin  specimen,  the 
material  would  buckle  well  before).  The  important  aspect  to  notice  is  that  more  necking  and  less  bulging 
occur  in  the  nonlinear  than  the  linear  model  in  tension,  and  compression,  respectively.  With  respect  to  the 
difference  in  the  linear  model  among  Fig.  2a  and  2b,  2b  is  the  reverse  of  2a  (this  is  a  characteristic  of 
linear  theory).  However,  the  lack  of  this  symmetry  for  the  nonlinear  model  is  characteristic  of  the  Green- 
Lagrange  strain  tensor  and  is  caused  by  the  interplay  between  the  linear  and  quadratic  terms  in  the  tensor. 
In  this  paper,  the  results  from  the  linear  and  nonlinear  model  will  be  compared. 


Fig  3.  Phantom  membrane  stretches  (a)  baseline,  (b)  1  cm  stretch,  (c)  2  cm  stretch. 


2.3.  Nonlinear  Model  Experiments 

In  previous  work,  a  phantom  of 
simulated  skin  was  generated  from  using  a 
polyurethane  membrane  25  cm  long,  15  cm 
wide,  and  approximately  2  mm  thick  containing 
two  types  of  materials.  The  bulk  material,  representing  “normal”  skin,  was  chosen  to  be  Evergreen  10 
(Smooth-On,  Easton,  PA),  while  the  stiffer  Evergreen  50  was  used  to  create  a  5  cm  circular  inclusion  of 
full  depth  and  was  placed  at  the  center  of  the  phantom.  The  visible  surface  of  the  phantom  was  modified 
by  drawing  a  regular  grid  of  horizontal  and  vertical  lines  spaced  about  1  cm  apart  in  either  direction.  Fig. 
3  shows  the  skin  phantom  as  acquired  by  an  optical  camera  at  the  baseline,  1  cm,  and  2cm  stretch  levels. 
Separate  mechanical  tests  on  the  membrane  were  conducted  using  the  EnduraTEC  ELF  3200  material 
tester  (EnduraTEC  Systems  Group,  Bose  Corporation,  Minnetonka,  MN) 25 .  Table  1  reports  the  expected 
Young’s  Modulus  contrast  ratio  based  on  a  Hookean  solid  fitted  to  each  respective  stretch/strain  level. 
The  scale  of  lesion  to  field-of-view  size  is  the  anticipated  aspect  ratio  for  a  dermoscopic  system.  This 


Stretch 

CR  (Linear) 

CR  (Nonlinear) 

1  (0.5  cm) 

5.7 

5.0 

2  (1.0  cm) 

5.0 

4.8 

3  (1.5  cm) 

4.6 

4.7 

4  (2.0  cm) 

4.2 

4.4 

Table  1.  Expected  elastic  contrast  ratio. 
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image  data  was  used  as  an  input  to  the  MIE 
framework,  and  reconstructions  comparing  the 
linear  and  nonlinear  model  were  generated. 

2.4.  Three-Dimensional  Model  Experiments 

In  an  effort  to  test  the  MIE  algorithm  with 
more  realistic  images  of  the  skin  as  would  be 
expected  in  the  clinic,  a  simulation  study  was 
performed  on  an  image  obtained  from  the 
Dermatology  Image  Atlas  project 
(www.dermatlas.org,  “melanoma_l_0405 1 0”, 
contributed  by  Eric  Ehrsam,  M.D.,  Fig.  4a). 

Dermoscopic  images  provide  the  challenge  of 
relying  on  the  natural  pattern  instead  of  the  structured  grid  used  in  the  membrane  phantom  experiments 
(although  the  skin  could  be  printed  upon  with  an  ink  grid).  Previous  simulation  work  on  this  image  can 
be  found  in 25 . 

One  possible  critique  of  this  dermoscopic  framework  is  the  fear  that  underlying  layers  would 
confound  the  approach.  As  a  result,  we  have  modified  the  simulation  study  reported  in  25  (e.g.  Fig.  4b)  to 
be  more  realistic.  Fig.  5a-5c  shows  the  setup  of  the  new  simulation  study.  In  this  study,  six  10  cm  x  10 
cm  domains  were  constructed  which  had  different  layered  structures.  Three  domains  were  6mm  in  total 
thickness  (-4mm  epidermis/dermis,  2 mm  subcutaneous)  and  had  1cm  central  inclusions  varying  in  depth 
1mm,  2mm,  and  3mm,  respectively.  The  remaining  three  domains  were  the  same  except  that  the 
subcutaneous  layer  thickness  was  increased  to  7mm.  The  selection  of  subcutaneous  layer  thickness  32, 33 
and  stiffness  values34  were  based  on  the  literature.  Using  this  domain,  a  mechanical  aperture  that 
stretched  the  skin  was  simulated,  3D  displacements  were  calculated  (Fig.  5b),  and  then  projected  on  to  the 


Fig.  4.  (a)  Melanoma  lesion,  with  (b)  sample 

reconstruction. 


b 


mm 


D=1cm 


T  ▼ 


Dermis 


HJh=1,  2,  3  mm 


Subcutaneous  layer 


▼T 


H=4mm 


L=2,  7  mm 


Fig.  5.  (a)  Mesh  with  refinement  in  lesion  area,  (b)  simulated  compression  of  skin  by  device,  and  (c)  depth 
variation  (transect  T  is  shown  in  (a)). 


original  2D  mesh.  A  new  set  of  simulated  images  (deforming  Fig.  4)  was  generated  and  then  used  as 
input  to  the  2D  MIE  algorithm. 

2.5  Noninvasive  In  Vivo  Experiments 


In  addition  to  the  simulation  experiments,  a  second  camera  and  deformation  system  was 
constructed  using  a  Sony  XCD-X710CR  CCD  camera  with  a  Schneider  12  mm  lens  and  a  circular 
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polarizing  filter  for  the  optics  mounted  over  a  crude  spring-loaded  standard  set  of  pincers  that  were 
pressed  against  the  skin  surface  and  bound  with  commercial  adhesive.  Digital  images  (800  x  600 
resolution)  of  a  common  nevus  of  palpable  stiffness  that  was  2-3  mm  in  diameter  and  located  on  the  volar 
aspect  of  the  left  forearm  of  a  male  volunteer  were  acquired  in  both  relaxed  and  stretched  states.  It  was 
determined  that  the  reflection  and  scattering  of  ambient  fluorescent  lighting  interfered  with  this  particular 
setup  and  affected  the  input  image  quality,  so  an  artificial  grid  was  imposed  to  counteract  specularity. 


Stretch  1 


4 


3 


Stretch  2 


3 


Fig  6.  (a)  Linear  Hookean  reconstmctions  for  each  stretch  state,  and  (b)  the  geometrically  nonlinear  Hookean 
reconstmctions.  Contour  shows  inclusion  border. 


Images  were  acquired  in  the  baseline  and 
stretched  state  and  given  to  the  MIE  algorithm. 

3.  RESULTS 

Fig.  6  illustrates  the  reconstructions  of 
the  phantom  membrane  system  using  the  linear 
and  nonlinear  models,  respectively.  Fig.  7  reports 
how  the  objective  function  varies  between  the 
linear  and  nonlinear  model  reconstructions  of  Fig. 
6.  Fig.  8  demonstrates  the  dependence  of 
resolving  stiffness  (3:1,  6:1,  12:1)  on  depth  (1,2, 

Fig  7.  Objective  function  difference  between  linear  and  ^  mm)  f°r  the  18  combinations  (combinations  of 
nonlinear  models  over  all  stretches  and  resolutions.  ^  depths,  3  contrast  ratios,  and  2  subcutaneous 

layer  thicknesses).  Table  2  reports  an 
approximate  contrast  ratio  of  the  lesion-to-bulk  material  for  each  of  the  cases  shown  in  Fig.  8.  Fig  9. 
illustrates  the  results  from  the  in  vivo  experiment  conducted.  Fig.  9  shows  (a)  the  nevus  ,  (b)  the  finite 
element  grid,  and  (c,d)  reconstructed  elasticity  images.  Fig.  9c-d  illustrate  the  effect  of  incorporating 
increasing  degrees  of  a  priori  knowledge  of  the  actual  location  and  elasticity  distribution  into  the 
algorithm.  Fig.  9c  is  a  general  elasticity  reconstruction  of  the  nevus  with  lesion-related  initial  guess,  and 
Fig.  9d  maintains  the  spatial  prior  of  Fig.  9c  for  the  entire  duration  of  the  reconstruction  process. 


Strain  (%>  S  #  of  regions 


4.  DISCUSSION 


With  respect  to  Fig.  6,  the  most  important  features  to  note  are  a  very  subtle  improvement  in 
satisfying  the  inclusion  contour,  and  the  decrease  in  variability  of  the  contrast  in  the  geometric  nonlinear 
reconstructions.  This  is  consistent  with  the  behavior  in  Table  1  and  indicates  that  the  geometric  nonlinear 
reconstruction  is  performing  as  predicted.  One  thing  to  note  is  that  if  the  membrane  were  a  true  Hookean 
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2  mm  subcutaneous 

1  mm  2  mm 


7  mm  subcutaneous 

1  mm  2  mm 


3:1 


6:1 


OOOH 
OOP  000 
000 


Fig.  8.  3-D  effects  on  elasticity  images  with  varying  lesion  depth  and  subcutaneous  layer  depth.  Ratio  on  left 
shows  actual  contrast  while  colorbar  shows  reconstructed.  Fig.  4  shows  the  lesion  contour. 


2mm 

SQ 

1  mm 

2  mm 

3  mm 

3:1 

2.12 

2.55 

2.79 

6:1 

3.08 

3.76 

4.00 

12:1 

4.11 

4.71 

5.15 

7  mm 

SQ 

1  mm 

2  mm 

3  mm 

3:1 

1.93 

2.27 

2.34 

6:1 

2.87 

3.49 

3.55 

12:1 

4.03 

4.28 

4.53 

Table  2.  Lesion-to-bulk  contrast  ratio  at 
each  lesion  depth  (1mm,  2mm,  3  mm) 
and  target  contrast  level  (3:1,  6:1,  12:1) 
for  2mm  (top)  and  7mm  (bottom) 
subcutaneous  thicknesses. 

nonlinear  solid,  the  CR  would  be  the 
same  across  all  stretches,  i.e.  strain 
levels.  While  the  variability  across 
levels  was  reduced  between  the 
linear  and  nonlinear  models  for 
strain,  the  data  supports  a  more 
complex  constitutive  model  for  this 
phantom  system.  Fig.  7  also  demonstrates  some  consistency  with  the  extension  to  a  nonlinear  model. 
Here  we  see  that  at  the  largest  stretch  values,  the  difference  between  the  two  models  is  highest  with  the 
nonlinear  model  outperforming  the  linear. 

With  respect  to  Fig.  8,  and  Table  2,  despite  not  reaching  the  expected  contrast  ratio  values,  each  result 
localized  the  lesion  and  demonstrated  sensitivity  to  depth  and  contrast.  With  respect  to  contrast  ratio,  this 


Fig.  9.  (a)  Dermoscopic  scale  interrogation  (magenta  border  is  about 
3  mm  in  diameter),  (b)  reconstmction  mesh,  (c)  general  elasticity 
reconstruction,  (d)  incorporation  of  lesion  border  as  a  priori 
information. 
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simulation  was  more  challenging  than  previous  work.  The  results  indicate  that  varying  subcutaneous 
thicknesses  has  a  very  modest  effect  in  confounding  lesion  reconstructions  for  the  contrast  regime  shown. 
Fig.  8  does  suggest  that  at  some  contrast  levels  there  may  be  overlap  (e.g.  a  3mm,  3:1  lesion  with  2mm 
subcutaneous  layer  may  be  confused  with  a  1mm,  6:1  lesion  with  7mm  subcutaneous  layer).  This  could 
be  potentially  confounding  if  left  in  the  absence  of  other  information.  However,  Horejsi  et  al.  and  Moller 
et  al.  32, 33  do  provide  guidance  for  estimating  the  subcutaneous  tissue  thickness  in  a  general  population. 
Additionally,  they  use  a  relatively  inexpensive  optical  device  to  make  these  measurements.  There  is  little 
doubt  that  using  this  a  priori  information  would  reduce  variability  in  interpretation.  It  is  interesting  to 
speculate  that  if  a  particular  lesion  type  was  of  a  known  stiffness  a  priori  then  perhaps  the  reconstruction 
contrast  could  be  used  to  estimate  lesion  depth.  Looking  across  the  grids  of  a  particular  property  contrast, 
there  is  a  definitive  change  in  the  reconstruction  as  a  function  of  depth. 

Fig.  9  demonstrates  a  preliminary  attempt  of  applying  the  framework  to  an  in  vivo  system. 
Considering  there  was  very  little  control  over  this  system  and  no  ground  truth  was  available,  these  results 
qualitatively  confirm  the  potential  utility  of  MIE  in  evaluating  an  area  of  skin  with  a  region  of  differing 
elasticity.  While  the  results  in  Fig.  9  c,d  are  not  completely  satisfying  at  this  time,  it  should  be 
emphasized  that  the  camera-system  employed  a  coarse  resolution,  and  the  skin-to-device  coupling  system 
was  extremely  crude.  Given  this  somewhat  ad-hoc  experiment,  the  detection  of  any  anomaly  in  the 
correct  region  is  encouraging  for  this  approach.  It  was  also  interesting  to  note  the  increase  in  performance 
by  adding  the  lesion  extents  to  the  information  provided  to  the  algorithm. 

5.  CONCLUSIONS 

The  results  from  the  experiments  above  demonstrate  an  ability  to  use  the  MIE  framework  within 
the  dermoscopy  setting.  Methodological  concerns  regarding  the  use  of  geometric  nonlinear  models, 
three-dimensional  effects,  and  in  vivo  conditions  are  addressed.  The  results  indicate  that  while  geometric 
nonlinearities  do  modestly  affect  reconstructions,  nonlinear  material  models  may  be  needed  to  correct  for 
remaining  discrepancies.  A  modest  improvement  was  shown  using  the  geometric  nonlinear  model 
especially  at  large  strain  values  which  is  consistent  with  the  theoretical  development.  The  three- 
dimensional  effects  of  lesion  depth  and  varying  subcutaneous  layer  thicknesses  are  assessed.  Lesion 
depth  does  affect  contrast  ratio  whereas  subcutaneous  layers  affect  the  reconstructions  to  a  significantly 
lesser  degree.  Avenues  to  detect  lesion  depth  in  the  presence  of  a  prior  knowledge  of  a  lesion’s  stiffness 
may  have  surgical  implications.  Preliminary  in  vivo  work  suggests  that  lesion  characterization  is  possible 
although  specificity  and  sensitivity  of  the  method  await  further  study  and  will  need  a  considerably  more 
robust  acquisition  system. 
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ABSTRACT 

This  work  explores  an  inverse  problem  technique  of  extracting  soft  tissue  elasticity  information  via  nonrigid  model- 
based  image  registration.  The  algorithm  uses  the  elastic  properties  of  the  tissue  in  a  biomechanical  model  to  achieve 
maximal  similarity  between  image  data  acquired  under  different  states  of  loading.  A  framework  capable  of  handling 
fully  three-dimensional  models  and  image  data  has  been  recently  developed  utilizing  parallel  computing  and  iterative 
sparse  matrix  solvers.  For  this  preliminary  investigation,  a  series  of  simulation  experiments  with  clinical  image  data  of 
human  breast  are  used  to  test  the  robustness  of  the  algorithm  to  expected  mis-estimation  of  displacement  boundary 
conditions  encountered  in  real-world  situations.  Three  methods  of  automated  point  correspondence  are  also  examined  as 
means  of  generating  boundary  conditions  for  the  algorithm. 

Keywords:  elastography,  computational  modeling,  inverse  problem,  non-rigid  registration 


1.  INTRODUCTION 

The  characterization  of  the  mechanical  properties  of  tissue  is  an  important  potential  source  of  information  for  detection 
and  diagnosis  of  disease  processes.  For  example,  there  is  a  long-standing  clinical  appreciation  of  evaluating  tissue 
elasticity  through  palpation  in  the  physical  examination  and  correlating  differences  in  stiffness  with  possible 
pathological  states.  A  minimally  invasive  methodology  for  analyzing  tissue  deformation  through  imaging  and/or  image 
processing  techniques  is  a  central  goal  of  the  field  of  elastography  [1,11].  Application  of  such  methods  to  the 
interrogation  of  the  breast  [2,3],  skin  [4-6],  prostate  [7],  and  other  accessible  organ  systems  is  an  active  area  of  research. 

Many  of  the  current  elastography  methods  are  founded  in  ultrasound  (US)  and  magnetic  resonance  (MR)  imaging  and 
involve  the  estimation  of  induced  displacements  within  the  tissue  of  interest  to  infer  the  elasticity  distribution.  We  have 
recast  the  problem  as  a  physically-constrained  non-rigid  image  registration  utilizing  quasi-static  deformation  and  image 
similarity  metrics  that  reconstruct  the  spatial  distribution  of  elasticity  parameters.  This  technique  has  been  termed 
’modality-independent  elastography’  (MIE)  [8-10]  because  of  its  ability  to  handle  native  anatomical  images  from 
different  sources  with  relatively  simple  modifications  to  the  acquisition  procedure.  To  date,  data  from  MR,  X-ray 
computed  tomography  (CT),  and  digital  photography  have  been  used  to  drive  the  algorithm.  In  addition  to  the  necessary 
preparation  and  effort  involved  in  gathering  images,  the  other  major  input  to  this  reconstruction  method  is  the 
delineation  of  boundary  conditions  on  the  region  of  interest.  Because  this  process  currently  involves  varying  levels  of 
manual  interaction,  there  is  a  need  to  develop  a  protocol  that  is  both  effective  and  mostly  automated  for  determining 
point  correspondences.  The  objectives  of  this  work  are  to  test  the  effects  of  degradation  in  input  data  quality  on  the  end 
reconstruction  and  examine  candidate  methods  for  generation  of  displacement  boundary  conditions.  This  is  done  in  the 
context  of  evaluating  the  robustness  of  a  newly  realized  three-dimensional  version  of  MIE  by  performing  simulation 
experiments  with  randomized  noise  processes  and  comparing  the  fidelity  of  reconstructions  resulting  from  boundary 
conditions  generated  by  three  different  techniques  of  determining  surface  point  correspondence. 

2.  METHODS 


2.1  Elastographic  reconstruction  framework 

The  conceptual  framework  for  our  elastographic  reconstruction  has  been  previously  described  in  [6,8-10].  In  brief,  an 
image  of  a  tissue  of  interest  {source)  is  deformed  by  a  biomechanical  computer  model  and  compared  against  an  acquired 
image  of  the  same  tissue  in  a  mechanically  loaded  state  {target).  Iterative  updates  of  elasticity  parameters  to  the  model 
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are  performed  until  a  suitable  match  in  intramodal  image  similarity  is  achieved  in  a  least  squares  manner  to  satisfy  a 
non-linear  optimization  scheme.  This  process  can  be  classified  as  an  inverse  problem,  with  model-based  deformation  of 
the  source  image  representing  the  forward  problem.  The  three  major  components  of  the  algorithm  are  the  model,  image 
comparison,  and  optimization,  each  of  which  is  described  in  more  detail  below. 

The  partial  differential  equation  that  expresses  a  state  of  mechanical  equilibrium  can  be  written  as  [12]: 

V  -cr  =  0  (1) 


where  <ris  the  Cartesian  stress  tensor.  We  have  elected  to  model  the  constitutive  tissue  behavior  using  Hooke’s  Law  of 
linear  elasticity,  which  states  that  the  strain  is  proportional  to  the  applied  stress,  and  assume  that  materials  are  isotropic 
and  incompressible  in  nature.  This  leads  to  the  formulation  of  Cauchy’s  Law  as 
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which  describes  the  constitutive  relationship  between  stress  and  strain  in  terms  of  the  elasticity  parameters  E  (Young’s 

£ 

modulus)  and  v  (Poisson’s  ratio).  The  shear  modulus  G  is  defined  as - . 

2(1  + v) 


A  finite  element  (FE)  representation  of  the  model  is  constructed  from  the  source  image  and  assigned  appropriate 
boundary  conditions  based  on  estimated  displacement  or  stress.  The  standard  Galerkin  method  of  weighted  residuals 
[13]  is  used  to  construct  and  solve  the  system,  which  yields  a  set  of  displacements  that  are  used  to  deform  the  source 
image.  This  model-deformed  image  is  then  compared  to  the  target  using  an  intensity-based  image  similarity  calculated 
for  a  series  of  voxel  groupings  determined  by  a  downsampling  of  the  image  set  overlap.  The  correlation  coefficient 
(CC)  [14]  is  the  method  of  choice,  as  it  has  empirically  demonstrated  superior  performance  over  other  metrics  such  as 
the  sum  of  squared  differences  and  normalized  mutual  information. 

The  elasticity  parameter  optimization  can  be  written  as  the  minimization  of  the  least  squares  error  objective  function 


Y  =  \s. 


TRUE 


f  EST 


(3) 


where  STrue  is  the  set  of  similarity  values  achieved  when  comparing  the  target  image  to  itself,  SEst  is  the  similarity 
between  the  model-deformed  source  and  the  target  images  using  current  estimates  of  the  elastic  modulus  distribution, 
and  |*|  denotes  the  vector  L2  norm.  Note  that  by  definition,  STRUE  for  CC  has  a  constant  value  of  1.  Using  a  Levenberg- 
Marquardt  approach,  the  residual  form  of  equation  (3)  becomes 


[jr/  +  a/]!A£}=[jr|sr„l,I  -Slsr} 


(4) 
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where  J  =  5SEst/^E  is  the  Jacobian  matrix  and  the  regularization  parameter  a  is  determined  using  the  methods  described 
in  [15].  Modulus  values  are  updated  by  AE  until  an  error  tolerance  is  reached  or  a  maximum  number  of  iterations  have 
been  completed.  Spatial  averaging  of  elasticity  values  in  the  model  and  solution  relaxation  between  iterations  are  also 
utilized  to  improve  the  stability  of  the  optimization. 

It  should  be  noted  that  the  size  of  the  Jacobian  matrix  is  dependent  on  the  number  of  material  properties  to  be 
reconstructed,  with  each  column  requiring  a  forward  solve  of  the  FE  model.  For  the  general  lesion  detection  problem,  a 
fine  discretization  of  the  mesh  requires  many  solutions  such  that  the  building  of  this  matrix  consumes  a  considerable 
portion  of  computational  resources.  This  fact  is  exacerbated  with  the  use  of  three-dimensional  data  and  necessitates  a 
parallelized  system.  Recent  work  using  the  Portable  Extensible  Toolkit  for  Scientific  Computation  (PETSc)  toolkit 
[16,17]  has  provided  the  necessary  coding  base  for  interfacing  sparse  matrix  system  solvers  with  a  C/C++  optimization 
framework  in  order  to  supersede  the  original  MATLAB/FORTRAN  framework.  This  new  version  is  designed  to  scale 
readily  between  the  complexity  of  the  model  and  the  number  of  processors  available;  it  has  been  tested  on  a 
homogeneous  cluster  of  ten  processors,  with  further  active  development  taking  place  in  conjunction  with  the  Vanderbilt 
ACCRE  project,  which  houses  hundreds  of  CPUs. 

2.2  Simulation  experiment  setup 

A  CT  volume  of  a  human  breast,  obtained  from  UC-Davis  Dept,  of  Radiology,  was  used  as  the  source  image  (256  x  256 
x  130,  0.6mm  x  0.6mm  x  0.6mm  voxel  spacing)  for  the  remainder  of  this  work.  The  surface  of  the  breast  was  segmented 
(ANALYZE  6.0,  Mayo  Cline,  Rochester,  MN)  to  create  a  three-dimensional  mesh  composed  of  39,013  nodes  connected 
as  214,163  tetrahedral  elements.  In  order  to  ensure  initial  data  fidelity  for  reconstruction  experiments,  an  idealized  target 
image  volume  and  gold  standard  boundary  condition  set  were  created.  A  2-cm  spherical  tumor  was  implanted  in  the 
center  of  the  mesh  by  assigning  a  stiff  modulus  value  to  the  member  elements  that  was  six  times  higher  than  the 
surrounding  material  [18].  Tissue  deformation  from  the  inflation  of  a  rectangular  air  bladder  against  the  lateral  surface  of 
the  breast  was  numerically  simulated  to  qualitatively  match  observed  mechanically  loading  of  an  existing  device  on  a 
breast-mimicking  phantom  of  polyvinyl  alcohol  cryogel.  The  stress  distribution  over  a  rectangular  contact  area  was 
modeled  as  the  cross-section  of  a  Gaussian  pressure  field  with  its  maximum  value  located  at  the  center  of  the  bladder; 
the  base  of  the  breast  was  fixed  in  place  as  if  pinned  to  the  chest  wall.  The  deformation  field  throughout  the  domain  was 
calculated  using  a  direct  forward  solve  of  the  finite  element  model  and  then  applied  to  the  intensity  field  of  the  source 
image  to  create  a  target  volume.  Displacements  at  the  surface  nodes  were  used  to  make  a  final  description  comprised  of 
all  Type  I  (Dirichlet)  boundary  conditions.  Figure  1  below  illustrates  the  setup  of  the  data  used. 

All  reconstructions  were  performed  using  a  priori  knowledge  of  the  location  and  size  of  the  inclusion  in  order  to  limit 
the  scope  of  the  problem  (e.g.  the  expense  of  the  Jacobian  matrix)  to  a  two-material  discrimination  of  relative  stiffness 
(elastic  contrast).  Having  a  defined  physical  model  and  synthetic  image  comparison  also  allows  for  examination  of  the 
optimization  behavior  separately  from  the  other  MIE  components  in  order  to  best  evaluate  the  effect  of  input 
inaccuracies  on  the  final  elasticity  distribution.  The  reconstruction  algorithm  begins  by  assigning  a  homogeneous 
elasticity  distribution,  with  Poisson’s  ratio  held  constant  at  v  =  0.485  to  represent  a  nearly  incompressible  material.  For 
this  data  set,  733  similarity  zones  were  discretized  from  the  target  image  volume. 


Figure  1.  Surface  renderings  of  CT  breast  volume  used  in  MIE  simulation  experiments.  From  left  to  right:  source  image, 
finite  element  mesh,  and  target  image  with  deformation  created  by  presumed  inflation  of  an  air  bladder. 
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2.3  Testing  robustness  of  the  algorithm 

The  current  method  of  selecting  boundary  conditions  as  derived  from  experiences  in  2D  work  requires  manual 
interaction  to  guide  or  correct  point  correspondence  for  every  surface  node.  Assuming  that  visible  markers  are  available 
in  an  image,  but  that  an  input  device  (e.g.  a  mouse)  is  needed  to  identify  the  specific  coordinates,  this  introduces  an 
operator-dependent  noise  process  in  localizing  any  given  point.  The  cumulative  effect  of  these  inaccuracies  is 
propagated  from  the  model  to  the  image  deformation  and  then  the  similarity  measurements.  For  a  given  reconstruction 
experiment,  the  gold  standard  boundary  condition  set  was  systematically  disrupted  by  adding  a  Gaussian  randomized 
vector  of  a  particular  length  (0.1,  0.2,  0.5,  1.0  or  2.0  voxel  units).  Figure  2  shows  an  example  of  the  distortion  caused  by 
the  applied  noise. 


Figure  2.  Example  of  distortion  due  to  additive  randomized  error.  The  gold  standard  boundary  conditions  used  to  generate 
a  controlled  deformation  produce  the  mesh  shown  on  the  left.  For  effect,  the  2.0  voxel  unit  noise  is  show  on  the  right. 


2.4  Testing  automated  boundary  condition  generating  methods 

For  this  work,  three  methods  of  surface  registration  and  point  correspondence  were  considered  for  a  more  automated 
method  of  determining  boundary  conditions  for  the  reconstruction  algorithm.  Two  are  derived  from  surface  matching  of 
potential  energy  distributions,  and  the  other  is  a  free-form  warping. 

If  the  flow  of  a  substance  over  undeformed  and  deformed  breast  surfaces  is  taken  to  be  a  conserved  process,  then 
correspondence  can  be  achieved  by  matching  the  same  energy  deposition  between  the  source  and  target,  that  is,  the 
equivalent  level  sets.  In  a  physical  sense,  Laplace’s  equation  can  be  used  to  describe  this  type  of  movement  analogously 
to  steady-state  heat  distribution: 

V2®  =  0  (5) 


where  O  would  represent  the  temperature  over  a  given  region.  Similarly,  the  diffusion  equation  describes  the  change  in 
concentration  or  density  of  a  material  over  time  on  a  region: 

5®  „  o -r. 

—  =  aV  ®  (6) 

dt 


where  a  is  the  diffusivity  constant. 

These  partial  differential  equations  were  used  to  calculate  an  energy  distribution  from  the  nipple  area  to  the  chest  wall 
over  the  surface  of  a  breast.  Isocontours  of  particular  energy  values  were  then  extracted  from  each  surface  to  form  a  set 
of  connected  points.  The  symmetric  closest  point  method  described  by  [19]  was  used  to  determine  a  displacement  field 
from  which  point  correspondence  at  boundary  nodes  could  be  interpolated. 
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The  third  method  involves  thin-plate  spline  interpolation  [20]  to  determine  point  correspondence.  This  was  done  to 
consider  a  widely-used  method  of  non-rigid  transform  that  can  take  advantage  of  fiducial  information  that  should  be 
present  in  future  real-world  data  acquisition.  The  use  of  physical  markers  to  track  breast  surface  displacement  during  a 
deformation  also  defines  a  set  of  control  points  would  allow  the  displacement  boundary  conditions  for  MIE  to  be  simply 
interpolated  from  the  local  warping.  In  these  simulation  experiments,  40  surface  nodes  of  known  correspondence  in 
each  of  the  image  volumes  (due  to  the  controlled  deformation)  served  as  fiducials. 

2.5  Evaluation  of  reconstructions 

Evaluation  of  the  reconstruction  results  is  performed  by  calculating  the  ratio  of  the  elasticity  of  the  inclusion  to  the  rest 
of  the  breast  for  the  distribution  that  yields  the  minimal  objective  function  value  over  the  course  of  optimization.  The 
robustness  of  the  MIE  algorithm  was  tested  with  four  trials  at  each  of  the  magnitudes  of  randomized  vectors  (described 
above  in  Section  2.3),  and  the  reconstruction  results  were  averaged  to  determine  a  trend  and  possible  threshold  of  noise 
tolerance  for  the  algorithm. 

For  the  automated  boundary  condition  generation  methods,  a  forward  mapping  of  the  objective  function  space  was 
calculated  to  determine  a  theoretical  optimum  to  the  reconstruction.  This  was  done  by  calculating  the  similarity  values 
for  model-based  image  deformations  created  by  adjusting  the  elasticity  contrast  of  the  inclusion  over  a  range  of  0.5:1  to 
30:1.  An  interpolating  curve  was  fit  and  the  minimum  objective  function  value  and  associated  elasticity  contrast  were 
extracted. 


3.  RESULTS 


3.1  Robustness  of  algorithm  to  boundary  condition  noise 

The  following  tables  summarize  the  effects  of  additive  noise  of  a  particular  magnitude  to  the  gold  standard  boundary 
condition  set.  As  the  magnitude  of  the  applied  randomized  vectors  increased,  a  dramatic  increase  in  the  minimum 
objective  function  value  is  observed.  Additionally,  changes  in  the  reconstructed  elasticity  contrast  indicate  that  a  cutoff 
exists  in  the  ability  of  the  algorithm  to  achieve  a  successful  result  (recall  that  the  known  correct  ratio  is  6:1)  for 
disruption  by  vectors  of  length  0.5  voxel  units  or  higher. 


Table  1.  Effect  of  applied  random  boundary  condition  noise  on  objective  function  space  and  reconstructed  elasticity 
contrast  ratio. 


Randomized 
vector  magnitude 
(voxel  units) 

Mean  optimal 
objective  function 
value 

Mean  optimal 
elasticity  contrast 
value 

0.1 

2.85  ±0.0382 

5.62  ±0.421 

0.2 

10.1  ±0.367 

5.70  ±0.588 

0.5 

60.1  ±4.19 

2.36  ±0.393 

1.0 

80.2  ±0.561 

2.47  ±  0.266 

2.0 

104  ±3.42 

2.17  ±0.422 
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3.2  Use  of  automated  point  correspondence 

Figure  3  depicts  the  deformation  fields  found  by  the  various  automated  methods  which  were  converted  into  a  boundary 
condition  sets  and  run  through  the  reconstruction  algorithm.  Qualitatively,  the  displacements  found  by  the  diffusion 
method  are  quite  different  from  the  gold  standard  set,  while  the  results  from  the  solution  of  Laplace’s  equation  and  the 
thin-plate  spline  interpolation  appear  to  be  more  satisfactory.  The  mean  target  registration  error  of  the  three  methods 
confirms  this  with  the  spline  having  the  best  performance  (0.26  mm),  the  Laplace  method  next  (1.0  mm),  and  the 
diffusion  method  being  the  worst  (2.0  mm).  Inspection  of  Figure  4  further  demonstrates  that  the  imposition  of  an  inexact 
boundary  condition  set  on  the  model  has  a  distinct  effect  on  the  optimization  by  shifting  the  minimum  value  to  a 
different  position.  A  comparison  of  the  fit  with  the  reconstruction  in  both  objective  function  value  and  elasticity  contrast 
is  provided  in  Table  3  below  and  indicates  that  the  algorithm  is  mostly  in  agreement  with  the  predicted  values  for  the 
Laplace  and  thin-plate  spline  methods  but  not  as  well  for  the  diffusion  method. 


Figure  3.  Three  candidate  automated  methods  for  MIE  boundary  condition  generation.  Top  row,  from  left  to  right:  surface 
deformations  calculated  from  diffusion  energy  matching,  Laplace  solution  energy,  and  thin-plate  spline  interpolation. 
Bottom  row:  target  registration  error  (TRE)  distribution  for  each  method  when  compared  against  the  gold  standard  of 
known  correspondence.  The  diffusion-based  mesh  is  both  qualitatively  and  quantitatively  the  worst  performer.  The 
Laplace  solution  appears  to  capture  the  shape  of  the  bladder  indentation  more  precisely,  but  the  thin-plate  spline  has  the 
best  overall  accuracy  in  determining  point  correspondence. 


cbjeclrve  funcUMi  valiM 
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Figure  4.  Mappings  of  objective  function  value  vs.  elasticity  contrast  ratio  (tumoribreast)  as  affected  by  the  boundary 
condition  sets  generated  by  different  automated  methods  of  surface  point  correspondence.  The  minimum  value  of  each 
curve  corresponds  to  the  optimal  elasticity  contrast  that  can  be  achieved  by  the  algorithm  when  constrained  by  the 
inaccuracies  of  the  methods:  (a)  diffusion,  (b)  Laplace,  and  (c)  thin-plate  spline  interpolations. 


Table  3.  Comparison  of  automated  point  correspondence  methods  on  MIE  reconstruction  quality.  Predicted  values  are 
found  from  the  minimum  point  of  the  curves  shown  in  Figure  4. 


Method 

Predicted  minimum 
objective  function 
value 

Reconstructed 
minimum  objective 
function  value 

Predicted 
optimal  elasticity 
contrast 

Reconstructed 
optimal  elasticity 
contrast 

Diffusion 

58.5 

58.8 

30.0 

12.6 

Laplace 

5.59 

5.59 

9.55 

10.3 

Thin-plate  spline 

12.3 

12.3 

5.55 

5.66 
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4.  DISCUSSION 

The  results  of  the  boundary  condition  noise  experiment  are  interesting  because  they  indicate  that  improper  localization 
of  boundary  points  greater  than  or  equal  to  0.5  units  of  voxel  spacing  can  introduce  significant  error  to  the  reconstruction 
process  and  impair  its  ability  to  characterize  the  underlying  elasticity  distribution.  This  is  a  similar  result  to  prior  work 
done  in  two-dimensional  systems  in  which  successful  reconstructions  correlated  to  boundary  condition  selection  error 
limited  to  half  a  pixel  length  [21].  It  also  confirms  that  randomizing  the  vectors  is  a  significant  challenge  to  the  algorithm 
because  it  introduces  highly  non-physical  deformations  that  cause  backlash  in  the  finite  element  mesh  and  other 
numerical  anomalies. 

The  implausibility  of  performing  manual  selection  on  all  6,319  boundary  nodes  underscores  the  importance  of  finding  an 
automated  method  for  determining  point  correspondences,  especially  at  less  than  0.5  voxel  units  of  error.  In  these 
simulation  experiments,  energy  matching  from  the  solutions  of  the  diffusion  and  Laplace  equations  yield  boundary 
condition  sets  that  are  inadequate  for  reconstructing  a  proper  elasticity  contrast.  This  can  be  partly  explained  because  the 
mean  errors  of  those  surface  registration  techniques  (as  compared  to  the  gold  standard)  are  approximately  3.3  and  1.7 
voxel  units,  respectively,  which  based  on  the  randomized  trials  were  magnitudes  too  large  for  the  algorithm  to  handle. 
The  diffusion-based  boundary  conditions  also  proved  more  difficult  to  obtain  a  stable  solution  for  in  the  model,  which 
probably  contributed  to  the  mismatch  in  reconstructed  elasticity  contrast.  However,  the  results  obtained  from 
reconstructions  using  the  thin-plate  spline  method  are  encouraging  because  the  mean  error  was  0.43  voxel  units.  The 
reconstruction  behavior  in  that  case  was  consistent  with  the  predicted  objective  function  space  and  the  optimal  elasticity 
contrast  was  found  to  be  within  6%  of  the  true  value.  This  preliminary  result  appears  to  identify  the  use  of  thin-plate 
spline  interpolation  as  a  strong  candidate  for  generating  boundary  conditions  for  MIE.  The  use  of  40  control  points  is 
also  seen  as  a  reasonable  choice  for  data  acquisition  and  processing  in  order  to  capture  the  extent  of  expected 
deformation  processes. 


5.  CONCLUSIONS 

In  this  work,  we  have  demonstrated  the  effects  of  inaccuracies  in  boundary  condition  determination  on  an  elastography 
method  that  maximizes  the  similarity  between  images  of  a  tissue  of  interest  acquired  under  two  different  states  of 
mechanical  loading.  In  order  to  characterize  the  robustness  of  the  current  version  of  this  method,  which  has  been 
updated  to  handle  three-dimensional  data  in  a  parallelized  scheme,  randomized  vectors  were  applied  to  distort  a  gold 
standard  boundary  condition  set.  The  results  were  used  to  determine  a  threshold  of  accuracy  needed  in  order  to  still 
achieve  an  accurate  reconstruction.  In  order  to  streamline  the  pre-processing  involved  in  the  algorithm,  three  methods  of 
automated  point  correspondence  were  evaluated.  The  success  of  these  methods  correlated  with  their  mean  error  (relative 
to  the  true  displacements)  meeting  the  putative  cutoff,  and  initial  results  indicate  that  established  techniques  such  as  thin- 
plate  splines  hold  promise  for  generating  boundary  conditions. 
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ABSTRACT 

Recent  advances  in  breast  cancer  imaging  have  generated  new  ways  to  characterize  the  disease.  Many  analysis 
techniques  require  a  method  for  determining  correspondence  between  a  pendant  breast  surface  before  and  after  a 
deformation.  In  this  paper,  an  automated  point  correspondence  method  that  uses  the  surface  Laplacian  or  the  diffusion 
equation  coupled  to  an  isocontour  matching  and  interpolation  scheme  are  presented.  This  method  is  compared  to  a  TPS 
interpolation  of  surface  displacements  tracked  by  fiducial  markers.  The  correspondence  methods  are  tested  on  two 
realistic  finite  element  simulations  of  a  breast  deformation  and  on  a  breast  phantom.  The  Laplace  correspondence 
method  resulted  in  a  mean  TRE  ranging  from  1.0  to  7.7  mm  for  deformations  ranging  from  13  to  33  mm,  outperforming 
the  diffusion  method.  The  TPS  method,  in  part  because  it  utilizes  fiducial  information,  performed  better  than  the 
Laplace  method,  with  mean  TRE  ranging  from  0.3  to  1.9  mm  for  the  same  range  of  deformations.  The  Laplace  and  TPS 
methods  have  the  potential  to  be  used  by  analyses  requiring  point  correspondence  between  deforming  surfaces. 

Keywords:  Registration,  non-rigid,  breast,  deformation,  correspondence,  surface  Laplacian,  diffusion,  finite  element, 
thin-plate  spline,  interpolation 


1.  INTRODUCTION 

As  breast  cancer  is  estimated  to  kill  over  40,000  women  and  be  diagnosed  in  more  than  178,000  in  2007  [1],  the 
detection  and  treatment  of  breast  cancer  is  an  important  area  of  scientific  research.  Many  novel  techniques  to  aid  in 
tumor  detection  are  being  developed  that  exploit  the  difference  in  physical  properties  between  healthy  and  cancerous 
tissue.  Some  of  these  techniques  measure  the  optical,  electrical,  or  elastic  properties  of  tissue,  including  near-infrared 
tomography  [6],  electrical  impendence  tomography  (EIT)  [7],  ultrasound  elastography  (USE)  [8],  magnetic  resonance 
elastography  (MRE)  [9],  and  in  particular,  modality- independent  elastography  (MIE)  [2,3]. 

MIE  is  a  reconstruction  algorithm  for  elasticity  imaging  that  can  be  used  for  detecting  breast  tumors.  It  involves 
imaging  a  pendent  breast  before  and  after  a  compression  and  using  these  images  to  reconstruct  the  elastic  properties  of 
the  tissue  using  a  nonlinear  optimization  framework,  computer  models  of  soft-tissue  deformation,  and  standard  measures 
of  image  similarity.  Unique  to  MIE  is  its  ability  to  utilize  images  from  any  modality  such  as  MRI  or  CT,  as  well  as  its 
usage  of  image  similarity  measures  that  make  direct  displacement  measurements  unnecessary. 

One  requirement  of  MIE  is  an  automated  method  of  finding  point  correspondence  between  the  pendent  breast  surfaces 
before  and  after  compression,  needed  to  specify  the  boundary  conditions  for  the  elasticity  model.  As  the  breast  is 
composed  of  soft  tissue  that  deforms  non-rigidly,  standard  rigid  registration  methods  cannot  be  applied.  Previous  work 
in  non-rigid  registration  includes  using  splines  and  FEM  models  [11],  as  well  as  point-based  methods  such  as  the 
symmetric  closest  point  (SCP)  algorithm  [10]. 

In  this  paper,  two  automated  methods  that  use  the  Laplace  and  diffusion  equations  to  establish  point  correspondence 
between  deformed  breast  surfaces  were  developed  and  compared  to  a  standard  thin-plate  spline  (TPS)  interpolation 
method  [4]. 


michael.i.miga@vanderbilt.edu 


53 


2.  METHODOLOGY 

2.1  Laplace  and  diffusion  methods  of  finding  point  correspondence 

A  major  investigative  task  of  this  work  was  to  evaluate  whether  the  energy  distributions  modeled  by  a  partial  differential 
equation  (PDE)  over  an  undeformed  {source)  surface  and  a  deformed  {target)  surface  can  be  used  to  find  the 
correspondence  between  the  two  surfaces.  In  this  method,  the  Laplace  and  diffusion  equations  were  independently 
solved  over  the  source  and  target  meshes  using  the  finite  element  method  (FEM).  Laplace’s  equation  is  most  commonly 
used  to  describe  potential  flow  problems  such  as  in  thermal,  fluid,  and  electrostatic  systems  and  is  given  by 

V£g7*:i-Q  (1) 

where  ®  represents  the  potential  and  a  describes  the  spatially  varying  conductivity.  The  diffusion  equation  which  allows 
a  time-varying  potential  is  given  by 


-  =  ¥-(<*?#}  (2) 
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where  ®  represents  the  potential  and  a  is  the  diffusion  coefficient.  Let  Osource  refer  to  the  solution  to  the  Laplace  or 
diffusion  equation  over  the  source  surface,  and  let  ®target  refer  to  the  solution  over  the  target.  The  basic  premise  is  that 
the  potential  field  distributed  over  the  source  and  target  surfaces  as  calculated  by  the  Laplace  or  diffusion  equation  will 
provide  information  about  the  correspondence  between  the  source  and  target  surfaces. 

To  solve  the  equations,  Dirichlet  boundary  conditions  were  set  to  simulate  potential  flow  from  the  nipple  area  to  the 
chest  wall  over  the  surface  of  a  pendent  breast  (specifically,  nodes  in  the  nipple  and  chestwall  area  were  given  boundary 
values  of  1  and  0,  respectively).  To  solve  equation  1,  a  Galerkin  finite  element  method  is  used  whereby  the  equations 
are  expressed  along  the  surface  orientation  (g=1).  To  solve  equation  2,  a  similar  scheme  was  used  for  handling  the 
spatial  component  of  the  PDE  and  a  fully  implicit  backwards  Euler  scheme  was  used  for  time-stepping.  In  the  case  of 
equation  2,  a  no-flux  condition  was  prescribed  at  the  chest  wall,  and  the  potential  field  was  allowed  to  propagate  from 
the  nipple  (a=l).  In  this  calculation,  time-stepping  was  stopped  once  the  potential  field  reached  the  chest  wall. 

After  the  Laplace  or  diffusion  equation  was  solved  over  the  source  and  target  surfaces,  the  solutions  were  used  to 
establish  correspondence  between  the  source  and  target  nodes.  This  involved  two  distinct  processes:  finding 
correspondence  between  isocontours  of  <Dsource  and  ®target  and  then  “interpolating”  that  correspondence  back  to  every 
source  node  on  the  mesh.  In  the  first  step,  isocontours  were  extracted  from  ®source  and  Otarget  for  a  set  of  selected 
isovalues.  The  correspondence  between  the  source  and  target  isocontours  was  determined  by  aligning  the  contours  by 
their  centroids  and  using  the  SCP  algorithm.  In  the  second  step,  the  displacement  vectors  at  the  source  isocontours  were 
interpolated  to  all  source  nodes  using  a  thin-plate  spline.  The  final  correspondence  was  found  by  adding  these 
displacements  to  the  source  nodes  to  get  the  location  of  the  corresponding  point  on  the  target  surface. 

The  method  can  be  summarized  in  the  following  steps: 

1.  Obtain  the  undeformed  source  mesh  and  deformed  target  mesh  that  define  a  breast  surface  before  and  after 
deformation. 

2.  Assign  boundary  conditions  at  nipple  and/or  chest  wall  nodes. 

3.  Solve  PDE  (diffusion  or  Laplace)  over  the  source  and  target  meshes  using  FEM. 

4.  Extract  isocontours  on  the  source  and  target  surfaces. 

5.  Determine  point  correspondence  between  source  and  target  isocontours  using  SCP. 

6.  Interpolate  displacements  at  source  isocontours  to  all  source  nodes. 

2.2  Using  thin-plate  spline  interpolation  to  find  point  correspondence 

One  advantage  of  the  PDE-based  correspondence  methods  is  that  they  do  not  explicitly  rely  on  external  markers  to 
constrain  the  matching  process.  However,  when  real-world  data  is  acquired,  fiducials  are  anticipated  to  be  available. 
Therefore,  TPS  interpolation  is  another  method  that  can  be  used  to  find  point  correspondence  [4].  Although  there  are 
many  different  methods  for  interpolation,  including  polynomial  splines  and  B-splines  [11],  TPS  interpolation  was  chosen 
in  part  because  it  does  not  require  a  regular  grid,  the  effects  of  changing  a  control  point  are  localized,  and  it  is  a  standard 
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method  that  has  been  successfully  used  in  many  non-rigid  registration  applications.  In  the  simulation  experiments 
described  below,  TPS  interpolation  was  used  to  find  point  correspondence  by  choosing  40  points  on  the  source  surface  to 
act  as  fiducials.  The  known  displacements  at  these  nodes  were  then  interpolated  to  all  surface  nodes  using  TPS. 

The  Laplace,  diffusion,  and  TPS  methods  for  finding  point  correspondence  described  above  were  tested  on  two 
simulation  data  sets  and  a  breast  phantom. 


2.3  Simulation  experiments 

To  perform  a  controlled  test  of  the  methods  described  above  on  a  breast-shaped  surface,  a  CT  image  volume  of  a 
pendant  breast  (courtesy  of  the  Dept,  of  Radiology,  University  of  Califomia-Davis)  was  segmented  to  create  a  source 
surface  consisting  of  6,313  points.  Two  types  of  deformations  were  simulated  by  assuming  different  contact  geometries 
of  an  air  bladder  being  inflated  against  the  breast  surface.  Circular  and  rectangular  cross-sectional  areas  of  a  Gaussian 
stress  distribution  positioned  along  the  lateral  aspect  of  the  breast  were  used  to  define  Type  2  boundary  conditions  for  a 
finite  element-based  deformation;  the  base  was  made  to  be  affixed  to  the  chest  wall.  The  displacement  solutions,  based 
on  a  three-dimensional  linear  elastic  model  of  a  Hookean  solid,  were  applied  to  create  the  target  surfaces  for  the  two 
simulations  (Figure  1). 


Figure  1.  Breast  surface  point  sets,  (a)  Source  surface  extracted  from  CT  volume  of  a  breast,  (b)  Target  surface 
generated  from  first  simulation  using  circular  cross-section  of  a  Gaussian  stress  distribution,  (c)  Target 
surface  generated  from  second  simulation  using  rectangular  cross-section  of  a  Gaussian  stress 
distribution.  The  correspondence  between  the  source  and  target  surfaces  was  determined  using  the 
Laplace,  diffusion,  and  TPS  methods. 

2.4  Phantom  experiments 

A  breast  phantom  was  constructed  to  test  the  point  correspondence  methods  with  real-world  data.  The  phantom  was 
fabricated  from  an  8%  w/v  solution  of  polyvinyl  alcohol  that  was  frozen  in  the  upper  half  of  a  2-liter  beverage  container 
for  16  hours.  After  8  hours  of  thawing,  thirty-four  1-mm  stainless  steel  ball  bearings  were  implanted  directly  under  the 
surface  of  the  resulting  cryogel  to  act  as  fiducials. 

The  phantom  was  then  imaged  inside  a  custom-built  rectangular  chamber  designed  to  deliver  compression  by  means  of 
an  air  bladder  placed  against  the  surface  of  the  phantom  (Figure  2).  CT  images  (512  x  512  x  174,  0.54  x  0.54  x  1  mm 
voxel  spacing)  were  acquired  with  the  phantom  at  three  different  states  of  mechanical  deformation  (undeformed,  50%  of 
maximum  bladder  pressure,  and  full  inflation).  Triangular  surface  meshes  were  obtained  by  semi-automatic 
segmentation  of  the  image  volumes  using  the  surface  extraction  tools  in  ANALYZE  6.0  (Mayo  Clinic,  Rochester,  MN), 
and  the  coordinates  of  the  fiducial  centroids  were  localized.  These  meshes  contained  approximately  8127,  6777,  and 
8260  nodes,  respectively. 


Figure  2.  Experimental  system  for  image  data  acquisition.  A  polyvinyl  alcohol  cryogel  is  placed  within  a 
Plexiglas  chamber  with  its  surfaces  held  in  place  against  the  walls.  Compression  is  delivered  through  an 
air  bladder  (arrow)  inflated  manually  through  a  bulb  adapted  from  a  standard  sphygmomanometer. 

The  Laplace,  diffusion,  and  TPS  methods  were  tested  on  the  phantom  surface  meshes.  For  the  TPS  method,  30  of  the 
fiducials  was  used  in  the  interpolation  and  the  four  remaining  fiducials  were  reserved  for  validation.  The  fiducials  used 
in  interpolation  and  validation  were  selected  such  that  the  distribution  for  both  groups  over  the  surface  was  roughly  even 
and  included  the  deformed  region. 

2.5  Validation 

In  order  to  assess  the  accuracy  of  the  simulation  and  phantom  experiments,  the  target  registration  error  (TRE)  was 
calculated.  The  TRE  measures  the  error  between  the  correspondence  determined  by  the  registration  method  and  the  true 
correspondence  [11].  For  the  simulation  experiments,  the  TRE  was  calculated  as  the  Euclidean  distance  between  the 
corresponding  target  points  determined  by  the  Laplace,  diffusion,  or  TPS  method  and  the  true  target  points.  Since  the 
true  correspondence  between  the  source  and  target  surfaces  was  known,  the  TRE  was  calculated  for  each  source  node, 
and  the  average  and  maximum  were  reported.  For  the  phantom  experiment,  the  TRE  was  calculated  using  the  centroids 
of  the  bead  fiducials  implanted  directly  under  the  phantom  surface.  Since  the  “gold  standard”  correspondence  was 
known  only  at  the  fiducials,  the  TRE  could  only  be  calculated  at  these  locations. 

In  addition,  since  one  crucial  step  in  both  the  Laplace  and  diffusion  methods  is  to  find  point  correspondence  between  the 
source  and  target  isocontours  (step  6  of  algorithm  summary),  we  evaluated  how  well  the  SCP  algorithm  performed  in 
this  step  for  the  simulation  data.  To  accomplish  this,  the  SCP  method  was  given  a  set  of  source  isocontours  and  their 
true  corresponding  target  contours,  and  the  error  (the  Euclidean  distance  between  the  true  target  point  and  the 
corresponding  target  point  determined  by  SCP)  was  calculated. 

3.  RESULTS 


3.1  Simulation  1  (circular  deformation  source) 

The  Laplace  and  diffusion  equations  were  solved  over  the  surfaces  generated  from  simulation  1  (cranial-caudal 
deformation  source  with  maximum  displacement  of  33  mm)  to  find  point  correspondence  between  the  source  and  target 
breast  surfaces.  For  comparison,  TPS  interpolation  using  40  simulated  fiducials  was  also  used  to  find  point 
correspondence.  The  accuracy  of  each  method  was  assessed  by  calculating  the  TRE  at  each  surface  node  (Figure  3). 
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Figure  3.  TRE  displayed  for  breast  simulation  1  (circular  deformation  source)  when  (a)  Laplace  equation, 

(b)  diffusion  equation,  and  (c)  TPS  interpolation  were  used  to  find  point  correspondence.  The  TPS 
method  resulted  in  the  lowest  error  overall  (mean  TRE  0.4  mm),  followed  by  the  Laplace  method 
(mean  TRE  2.3  mm)  and  diffusion  method  (max  TRE  4.5  mm).  The  highest  TRE  is  found  in  the 
deformed  region  when  the  Laplace  method  is  used  and  in  the  base  when  the  diffusion  method  is 
used. 

The  results  (Table  1)  indicated  that  the  Laplace  method  performed  more  accurately  overall  than  the  diffusion  method; 
however,  the  area  with  the  highest  amount  of  error  differed.  When  the  Laplace  method  was  used,  the  deformed  region 
had  the  highest  error,  whereas  when  the  diffusion  method  was  used,  the  area  farthest  from  the  diffusion  source  had  the 
highest  error.  (In  this  case,  since  the  diffusion  source  was  located  in  the  nipple  area,  the  highest  error  occurred  in  the 
chest  wall  region.)  The  TPS  interpolation  had  the  lowest  error  overall,  and  the  error  distribution  over  the  surface  was 
related  to  the  locations  of  the  simulated  fiducials. 

The  results  given  above  pertain  to  a  simulated  compression  with  a  maximum  displacement  of  approximately  33  mm. 
Since  this  amount  of  compression  may  be  larger  than  is  needed  for  many  applications  and  may  introduce  other  unwanted 
effects  for  MIE  due  to  non-linear  elastic  behavior,  the  point  correspondence  methods  were  also  tested  for  lesser  amounts 
of  compression.  The  TRE  for  different  amounts  of  compression  when  the  Laplace  method  was  used  to  find  point 
correspondence  is  shown  in  Figure  4.  The  TRE  appears  to  increase  linearly  with  respect  to  increasing  compression. 

The  mean  and  maximum  error  for  the  isocontour  point  correspondence  determined  by  the  SCP  algorithm  (detailed  in 
methods  section)  was  calculated  (Figure  5).  The  isocontour  correspondence  given  by  the  SCP  algorithm  had  a 
maximum  error  of  about  5  mm  for  the  maximum  compression  of  33  m. 


57 


Figure  2.  Maximum  TRE  (solid  line)  and  mean  TRE  (dotted  line)  when  the  Laplace  method  was  used  to  find  point 
correspondence  between  the  source  and  target  surfaces  at  different  levels  of  compression.  A  deformation 
of  100%  indicates  a  maximum  displacement  of  approximately  33  mm.  The  TRE  seems  to  increase 
linearly  with  respect  to  the  amount  of  deformation 


Fraction  of  deformation 

Figure  5.  Evaluation  of  the  accuracy  of  the  SCP  method  used  to  find  isocontour  correspondence.  The  mean  error  of 
the  SCP  method  (dotted  line)  is  compared  to  the  mean  error  of  the  Laplace  method  (solid  line)  when  SCP 
method  was  used  to  find  isoncontour  correspondence  for  simulation  1  data  at  different  levels  of 
compression.  (A  deformation  of  100%  indicates  a  maximum  displacement  of  approximately  33mm.) 
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3.2  Simulation  2  (rectangular  deformation  source) 

When  the  Laplace,  diffusion,  and  TPS  interpolation  methods  were  used  to  find  point  correspondence  between  the  breast 
surfaces  generated  by  simulation  2  (a  more  realistic  simulation  using  a  rectangular  deformation  source  with  a  maximum 
displacement  of  13  mm),  the  results  (Figure  6)  were  very  similar  to  those  from  simulation  1.  However,  the  TRE  for  all 
three  methods  (Table  1)  was  slightly  lower,  possibly  due  to  the  lower  degree  of  compression  simulated. 


meters 


(C) 

Figure  3.  TRE  displayed  for  breast  simulation  2  (rectangular  deformation  source)  when  (a)  Laplace  equation, 
(b)  diffusion  equation,  and  (c)  TPS  interpolation  were  used  to  find  point  correspondence.  The  TPS 
method  resulted  in  the  lowest  error  overall  (mean  TRE  0.3  mm),  followed  by  the  Laplace  method 
(mean  TRE  1 .0  mm)  and  diffusion  method  (mean  TRE  2.1  mm).  The  results  are  similar  to  that  of 
simulation  1 
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3.3  Phantom 

The  Laplace  and  diffusion  methods  were  used  to  determine  point  correspondence  between  the  noncompressed  and 
compressed  surfaces  of  a  breast  phantom.  The  results  were  validated  by  calculating  the  TRE  at  34  fiducials  located 
directly  below  the  surface  of  the  phantom.  For  comparison,  TPS  was  used  to  interpolate  the  displacements  of  30 
fiducials  to  all  surface  nodes,  and  the  TRE  was  calculated  using  the  4  remaining  fiducials. 

The  results  for  a  50  and  100%  compression  (with  a  maximum  displacements  of  about  20mm  and  36  mm,  respectively) 
are  shown  in  Table  2.  As  in  the  simulations,  the  Laplace  method  performed  better  overall  than  the  diffusion  method  and 
had  lower  TRE.  The  TRE  for  the  TPS  interpolation  was  lower  than  that  for  the  Laplace  and  diffusion  methods,  but 
varied  with  the  number  and  locations  of  fiducials  used  in  the  interpolation. 


Figure  4.  Breast  phantom  surfaces  (a)  before  compression,  (b)  at  50%  compression  with  maximum 
displacement  of  .020  m,  and  (c)  at  100%  compression  with  maximum  displacement  of  .036  m.)  Lines 
indicate  isocontours  at  different  values  of  k.  Black  nodes  at  the  nipple  and  base  indicate  the  nodes 
assigned  boundary  values. 
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Table  1.  TRE  for  different  point  correspondence  methods  tested  on  breast  surfaces  generated  from  simulation  1 
(point  deformation  source  with  max  displacement  of  33  mm)  and  simulation  2  (rectangular 
deformation  source  with  max  displacement  of  13  mm)  breast  surfaces.  The  Laplace  method 
outperformed  the  diffusion  method,  while  the  TPS  method  performed  best  of  all. 


Simulation  1 

Simulation  2 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Laplace 

14.6 

2.3 

4.6 

1.0 

Diffusion 

24.2 

4.5 

10.3 

2.1 

TPS  (40  fiducials) 

7.6 

0.4 

2.6 

0.3 

Table  2.  TRE  for  different  point  correspondence  methods  tested  on  breast  phantom  at  approximately  50%  and 
100%  compression,  with  max  displacements  of  20  and  36  mm,  respectively.  The  Laplace  method 
outperformed  the  diffusion  method,  while  the  TPS  method  performed  best  of  all. 


50%  Compression 

1 00%  Compression 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Laplace 

8.1 

3.5 

16.4 

7.7 

Diffusion 

11.9 

3.9 

19.4 

7.9 

TPS  * 

1.5 

1.1 

2.6 

1.9 

*  TPS  interpolation  using  30  fiducials;  4  fiducials  were  used  to  calculate  TRE. 


4.  DISCUSSION 

Of  the  three  registration  methods  evaluated,  the  TPS  method  consistently  outperformed  the  Laplace  and  diffusion 
methods  and  had  the  lowest  error  for  both  the  simulation  and  phantom  experiments.  However,  it  is  important  to  note 
that  a  comparison  of  the  PDE-based  methods  and  the  TPS  method  is  not  entirely  fair  since  the  TPS  method  relies  on 
fiducial  information  that  the  Laplace  and  diffusion  methods  do  not  require.  The  performance  of  the  TPS  method  is 
dependent  on  both  the  number  and  placement  of  these  fiducials.  These  results  indicate  that  30-  40  fiducials  with  an  even 
distribution  over  the  surface  should  be  sufficient  to  register  surfaces  (with  13-33  mm  deformations)  with  mean  TRE 
ranging  from  0.3  to  1.9  mm.  Although  further  studies  are  needed  to  determine  the  optimal  number  and  placement  of 
fiducials,  experience  suggests  that  increasing  the  number  of  fiducials  in  the  areas  with  greatest  deformation  increases 
registration  accuracy,  and  conversely,  lowering  the  number  of  fiducials  in  those  areas  causes  a  significant  decrease  in 
accuracy. 

The  results  indicate  that  the  Laplace  method  is  a  useable  surface  registration  method.  Although  the  Laplace  method  did 
not  perform  as  accurately  as  the  TPS  method,  it  has  the  advantage  of  not  requiring  fiducial  information.  However,  one 
of  the  challenges  of  the  Laplace  method  is  determining  the  regions  to  which  boundary  conditions  are  assigned.  Accurate 
selection  of  these  regions  is  important  because  the  implicit  correspondence  between  these  regions  is  used  by  the  Laplace 
equation  to  obtain  the  correspondence  for  the  rest  of  the  surface.  Lor  these  studies,  the  nipple  region  and  the  chest  wall 
boundary  regions  were  selected  manually,  further  studies  may  be  needed  to  find  a  method  to  automate  the  selection  of 
the  boundary  regions  and  to  evaluate  how  error  in  the  selection  of  these  regions  affects  the  final  registration  error. 

Although  the  diffusion  method  does  have  certain  advantages  over  the  Laplace  and  TPS  registration  methods,  several 
problems  prevent  it  from  becoming  a  viable  surface  registration  technique.  The  advantages  of  the  diffusion  method  are 
that  the  correspondence  near  to  the  diffusion  source  (in  this  case,  the  nipple)  is  relatively  accurate.  In  addition,  the 
diffusion  method  only  requires  boundary  conditions  to  be  set  in  one  region  (in  this  case,  the  nipple),  unlike  the  Laplace 
method,  which  requires  boundary  conditions  at  two  regions  (nipple  and  chest  wall  base),  and  the  TPS  method,  which 
requires  multiple  points  of  constraint  (at  34  fiducials). 
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However,  the  diffusion  method  does  not  appear  to  be  a  very  useable  surface  registration  method  for  the  following 
reasons:  the  results  indicate  a  substantial  amount  of  error,  the  registration  and  resulting  error  are  highly  dependent  on  the 
diffusion  parameters  chosen  (time  step  and  final  time  in  particular),  and  the  diffusion  parameters  must  be  manually 
adjusted  for  each  different  surface  mesh  since  there  is  no  automated  method  to  find  the  optimal  diffusion  parameters. 
Since  the  diffusion  described  by  the  PDE  is  by  definition  a  non  steady-state  process,  an  optimal  registration  requires  that 
the  diffusion  front  should  travel  over  the  entire  surface  between  the  nipple  and  base  and  stop  at  the  base  in  order  to 
assure  correspondence  for  as  much  as  the  surface  as  possible.  If  the  parameters  are  chosen  such  that  the  diffusion  front 
does  not  reach  the  base,  the  correspondence  for  the  regions  not  reached  by  the  diffusion  front  cannot  be  constrained  and 
must  be  interpolated  from  the  displacements  of  the  surrounding  regions.  Conversely,  if  the  diffusion  front  travels  for  too 
long  a  time,  the  solution  over  the  surface  approaches  saturation,  resulting  in  a  flat  gradient  and  lack  of  isocontours  from 
which  to  establish  correspondence.  Various  modifications  to  the  diffusion  method  employing  curvature  information  and 
using  different  diffusion  coefficients  were  tested,  but  none  was  very  successful.  Therefore,  the  sensitivity  of  the 
diffusion  method  to  parameters  and  substantial  amount  of  error  may  prevent  the  diffusion  method  from  being  a  viable 
surface  registration  method. 

The  TRE  measured  for  each  registration  technique  is  not  only  dependent  on  the  factors  described  above,  but  also  on  the 
amount  of  deformation  of  the  target  surface.  The  results  suggest  that  the  TRE  increases  roughly  linearly  with  the 
amount  of  deformation.  Using  the  simulation  and  phantom  data  presented  here,  one  may  be  able  to  estimate  the  range  of 
error  expected  when  one  of  the  described  methods  is  used  to  register  breast  surfaces  with  a  particular  amount  of 
compression.  Conversely,  the  maximum  amount  of  compression  that  will  yield  a  registration  within  a  given  error  bound 
can  be  roughly  estimated.  For  the  purposes  of  MIE,  realistic  compressions  will  be  in  the  range  of  1-2  cm. 

Another  factor  related  to  the  amount  of  compression  is  the  distribution  of  TRE  over  the  surface.  The  TRE  was  not 
evenly  distributed;  rather,  the  TRE  in  the  areas  of  greatest  deformation  tended  to  be  higher  than  the  TRE  elsewhere. 
Therefore,  the  mean  TRE  is  not  necessarily  the  best  measure  of  the  TRE  over  the  surface;  the  max  TRE  may  reflect  the 
error  in  the  deformed  regions  more  accurately. 

In  addition  to  the  evaluation  of  the  three  registration  methods,  the  performance  of  the  SCP  algorithm  was  evaluated  since 
the  matching  of  the  isocontours  extracted  from  the  source  and  target  surface  is  a  crucial  step  of  the  PDE-based 
registration  methods.  The  results  indicate  that  the  amount  of  error  the  SCP  algorithm  contributes  to  the  Laplace  and 
diffusion  methods  is  relatively  small  when  compared  to  the  total  TRE  (Figure  2). 

In  comparsion  to  previous  studies,  the  Laplace  method  outperformed  the  modified  SCP  method  implemented  by  Schuler, 
et  al.  The  data  generated  by  first  simulation  described  in  this  paper  was  also  used  to  test  the  modified  SCP  method,  and 
whereas  the  Laplace  method  had  a  maximum  error  of  14.6  mm  for  a  deformation  of  33  mm,  the  modified  SCP  method 
had  a  maximum  error  of  27.8  mm  [5]. 

MIE  is  one  application  that  may  use  the  registration  methods  described  in  this  paper,  in  this  case  to  determine  boundary 
conditions  for  its  elasticity  model.  Preliminary  studies  indicate  that  TPS  is  the  most  viable  registration  method,  the  error 
of  which  is  within  the  bounds  required  for  a  successful  elasticity  reconstruction  (approximately  0.3  mm).  The  mean 
error  for  the  Laplace  registration  method  exceeds  MIE’s  error  bounds,  and  although  the  target  boundary  conditions 
produced  by  the  Laplace  method  resulted  in  a  viable  mesh,  the  resulting  elasticity  reconstruction  contained  a 
considerable  amount  of  error.  The  diffusion  method  could  not  be  used  in  conjunction  with  MIE  because  of  the  extreme 
distortion  of  the  target  finite  element  mesh  generated  from  the  surface  registration. 

5.  CONCLUSION 

The  results  of  the  simulation  and  phantom  experiments  indicate  that  while  TPS  interpolation  is  the  most  accurate  surface 
registration  method  of  those  evaluated,  the  Laplace  method  may  be  a  viable  surface  registration  technique  if  fiducials  are 
not  available.  Although  the  TPS  method  consistently  outperformed  the  Laplace  method,  its  performance  is  dependent  on 
the  number  and  distribution  of  fiducials  available.  Both  the  Laplace  and  TPS  methods  have  been  used  in  MIE  to  register 
breast  surfaces  in  order  to  determine  boundary  conditions  for  its  elasticity  model.  In  addition  to  MIE,  the  Laplace  and 
TPS  methods  also  have  potential  to  be  used  for  non-rigid  registration  in  more  general  applications. 
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